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ABSTRACT 


Through-thickness  reinforcement  (stitching,  z-pins)  is  a  promising  technology  to  develop  fail-safe  load  bearing 
components  for  aeronautical  structures  and  lightweight  armor  and  combat  vehicles  with  superior  capabilities.  The 
through-thickness  reinforcement  enhances  the  performance  of  laminates  against  delamination  failure,  improving  impact 
and  damage  tolerance  and  rendering  stable  or  even  suppressing  crack  propagation.  The  through-thickness  reinforcement 
creates  large  zones  of  bridging  along  the  wake  of  the  cracks  developing  tractions  that  reduce  the  crack  driving  force. 

The  purpose  of  the  work  was  to  investigate  basic  aspects  of  the  mechanics  of  through-thickness  reinforced  laminates. 
Focus  has  been  on  two  problems. 

1)  The  definition  and  validation  of  approximate  mode  I  weight  functions  for  orthotropic  double  cantilever  beams.  The 
weight  functions  allow  the  problem  of  large  scale  bridging  delamination  of  composite  laminates  to  be  formulated  as 
integral  equations  overcoming  the  limitations  imposed  on  accuracy  by  beam  theory  approximations.  The  weight 
functions  have  been  deduced  from  the  corresponding  isotropic  result  using  an  orthotropy  rescaling  technique  and  have 
been  applied  to  investigate  the  influence  of  the  orthotropy  on  the  fracture  behavior  of  laminates  in  the  presence  of 
Dugdale’s  type  crack  bridging  mechanisms. 

2)  The  formulation  of  a  theoretical  model  to  investigate  the  effectiveness  of  a  through-thickness  reinforcement  in 
reducing  the  delamination-induced  degradation  of  the  dynamic  properties  of  delaminated  structures.  The  model  is  based 
on  the  first  order  shear  deformation  plate  theory  and  deals  with  plates  with  through-width  delaminations.  The  bridging 
mechanisms  are  modeled  as  a  uniform  distribution  of  linear  elastic  springs  that  oppose  the  relative  crack  displacements. 
The  model  has  been  applied  to  predict  the  free  vibrations  of  a  delaminated  cantilever  beam.  A  variety  of  solutions  and 
transitions  in  the  dynamic  response  has  been  found  on  varying  the  stiffness  of  the  ligaments.  An  application  to  a 
stitched  carbon-epoxy  laminate  for  the  aeronautical  industry  shows  that  low  area  fractions  of  through-thickness 
reinforcement  can  substantially  improve  the  dynamic  response  of  delaminated  structures  and  almost  entirely  restore 
their  integrity. 


KEYWORDS:  through-thickness  reinforcement;  advanced  composite  laminates,  delamination  fracture;  free 
vibrations;  weight  function  methods;  plate  theories. 
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TECHNICAL  REPORT 


BACKGROUND  AND  STATEMENT  OF  THE  PROBLEM 

The  purpose  of  the  work  was  to  investigate  basic  aspects  of  the  mechanics  of  through-thickness  reinforced  laminates. 
Through-thickness  reinforcement  (stitching,  z-pins)  is  a  promising  technology  to  develop  fail-safe  load  bearing 
components  for  aeronautical  structures  and  lightweight  armor  and  combat  vehicles  with  superior  capabilities.  A  wide 
body  of  theoretical  and  experimental  research  has  been  produced  over  the  last  ten  years  that  leaves  no  doubt  that  a 
through-thickness  reinforcement  substantially  enhances  the  performance  of  laminates  against  delamination  failure, 
improving  impact  and  damage  tolerance  and  rendering  stable  or  even  suppressing  crack  propagation  (see  references 
cited  in  [1-4]).  The  through-thickness  reinforcement  creates  large  zones  of  bridging  along  the  wake  of  the  cracks 
developing  tractions  that  oppose  the  relative  crack  displacement  and  reduce  the  driving  force  for  crack  propagation. 

In  this  project  focus  has  been  on  two  problems. 

1)  The  definition  and  validation  of approximate  mode  I  weight  functions  for  orthotropic  double  cantilever  beams. 

The  weight  functions  allow  the  problem  of  large  scale  bridging  delamination  of  composite  laminates  to  be  formulated 
as  integral  equations  overcoming  the  limitations  imposed  on  accuracy  by  beam  theory  approximations.  Such  an 
approach  will  permit  the  simulation  of  unusual  phenomena  of  crack  face  closure  and  crack  propagation  in  the  presence 
of  regions  of  contact  that  have  been  recently  observed  in  tests  performed  on  through-thickness  reinforced  laminates  and 
that  cannot  be  properly  described  by  approximate  theories. 

2)  The  formulation  of  a  theoretical  model  to  investigate  the  effectiveness  of  a  through-thickness  reinforcement  in 
reducing  the  delamination-induced  degradation  of  the  dynamic  properties  of  delaminated  structures. 

It  is  expected  that  a  through-thickness  reinforcement  would  reduce  the  degradation  of  the  dynamic  properties  of 
delaminated  structures  through  the  development  of  tractions  acting  along  the  surfaces  of  the  delamination  and  opposing 
the  relative  crack  displacements  (opening  and  sliding).  This  problem  has  been  investigated  under  this  contract  and  the 
effects  of  a  through-thickness  reinforcement  quantified  through  the  proposed  model. 


APPROACH 

1)  The  weight  functions  have  been  deduced  from  the  corresponding  isotropic  result  using  an  orthotropy  rescaling 
technique  and  numerical  results  to  match  asymptotic  solutions.  The  weight  functions  have  been  used  in  the  formulation 
of  a  nonlinear  fracture  mechanics  model  where  the  problem  has  been  posed  as  integral  equations  in  terms  of  crack  tip 
stress  intensity  factors,  (see  [1,3]  for  details) 

2)  A  model  based  on  the  first  order  shear  deformation  plate  theory  has  been  formulated  to  analyze  the  dynamic  response 
of  plates  in  the  presence  of  through-width  delaminations.  The  bridging  mechanisms  are  modeled  as  a  uniform 
distribution  of  linear  elastic  springs  that  oppose  the  relative  crack  displacements  and  develop  symmetric  restoring  forces 
during  vibration.  The  model  has  been  applied  to  predict  the  natural  frequencies  of  vibration  and  modal  shapes  of  a 
delaminated  cantilever  beam.  The  natural  frequencies  are  defined  through  the  solution  of  a  boundaty  eigenvalue 
problem.  An  analytical  solution  has  been  found  for  the  case  of  homogeneous  material,  negligible  shear  deformations, 
mid-plane  delamination  and  constrained  model  assumption.  A  numerical  procedure  (based  on  a  relaxation  technique) 
has  been  formulated  to  solve  general  cases  (see  [2,4]  for  details). 


RESULTS 

During  the  course  of  the  project  three  papers  have  been  published  (one  in  conference  proceedings  [4]  and  two  in 
international  journals  [1,2])  where  the  results  and  conclusions  of  the  work  have  been  presented  in  great  detail.  A  list  of 
the  main  results  obtained  under  this  contract  and  future  development  of  the  work  is  given  below: 

1 )  Mode  I  weight  functions  for  orthotropic  double  cantilever  beams. 

Mode  I  weight  functions  have  been  derived  and  validated  numerically  that  permit  a  single  model  to  treat  the 
progression  from  the  initiation  of  a  bridged  crack  to  its  propagation  as  a  long  crack  with  large  scale  bridging  in 
orthotropic  materials. 
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A  nonlinear  fracture  mechanics  model  has  been  formulated  to  analyze  crack  propagation  under  general  large  scale 
bridging  conditions  and  results  have  been  presented  for  the  case  of  bridging  mechanisms  characterized  by  a 
Dugdale’s  type  law. 

Based  on  the  results  obtained  under  this  contract,  work  is  currently  in  progress  on  the  definition  of  mode  II  weight 
functions  for  orthotropic  double  cantilever  beams.  The  weight  functions  will  allow  the  analysis  of  crack 
propagation  under  mode  I  (double  cantilever  beams),  mode  II  (end-notched  flexural  specimens)  and  mixed  mode 
(mixed-mode  bending  specimen)  conditions. 

The  procedure  proposed  in  [1]  will  be  extended  to  define  mode  I  weight  functions  for  dynamic  crack  propagation. 


2)  Dynamic  response  of  delaminated  beam-type  structures  in  the  presence  of  through-thickness  reinforcement 

A  theoretical  model  has  been  formulated  to  analyze  the  free  vibration  response  of  composite  laminate  beam-type 
structures  in  the  presence  of  linear-elastic  crack  bridging  mechanisms  acting  along  the  crack  surfaces. 

The  model  has  been  applied  to  the  case  of  a  cantilever  beam  and  a  variety  of  solutions  and  transitions  in  the 
dynamic  response  has  been  found  on  varying  the  stiffness  of  the  ligaments. 

An  application  to  a  stitched  carbon-epoxy  laminate  for  the  aeronautical  industry  has  shown  that  low  area  fractions 
of  through-thickness  reinforcement  can  substantially  improve  the  dynamic  response  of  delaminated  structures  and 
almost  entirely  restore  their  integrity;  in  particular  it  has  been  observed  that  the  opening  modes  are  completely 
avoided  in  the  presence  of  crack  bridging  mechanisms. 

Based  on  the  results  obtained  under  this  contract,  work  is  currently  in  progress  to  investigate  and  summarize  the 
effectiveness  of  a  through-thickness  reinforcement  for  different  boundary  conditions,  positions  and  lengths  of  the 
delamination  in  plates  vibrating  in  cylindrical  bending.  This  summary  will  be  used  to  guide  the  optimal  design  of 
through-thickness  reinforced  structures.  Future  development  will  be  in  the  investigation  of  the  dynamic  response  of 
delaminated  structures  in  the  presence  of  nonlinear  crack  bridging  mechanisms. 


DIVULGATION  OF  THE  RESULTS 

Apart  from  being  published  in  [1-4],  the  work  completed  during  the  project  has  been  presented  at  three  conferences: 

1 0th  Int.  Congress  of  Fracture,  ICF1 0 ,  December  2001 ,  Honolulu,  USA. 

15th  ASCE  Engineering  Mechanics  Conference ,  June  2002,  New  York,  USA 

14th  US  National  Congress  of  Theoretical  and  Applied  Mechanics ,  July  2002,  Blacksburg,  VA,  USA  (special 
symposium  on  Impact  of  Composites). 
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International  Journal  of  Engineering  Science,  2002,  in  press. 

Mode  I  Weight  Functions  for  an  Orthotropic  Double  Cantilever  Beam 

R.  Massabo  L.  Brandinelli ',  and  B.  N.  Cox2 
1  Department  of  Structural  and  Geotechnical  Engineering, 

University  of  Genova,  Genova,  16145,  Italy 
2  Rockwell  Scientific  Co.,  Thousand  Oaks,  CA  91360,  USA 

Abstract 

Approximate  weight  functions  are  proposed  and  validated  numerically  for  an  orthotropic  double 
cantilever  beam  loaded  in  mode  I.  They  define  the  stress  intensity  factor  at  the  crack  tip  due  to  a  pair 
of  point  forces  acting  on  the  crack  surfaces  and  have  been  deduced  from  the  corresponding  isotropic 
result  using  an  orthotropy  rescaling  technique.  The  weight  functions  allow  mode  I  large  scale  bridging 
problems  in  beams  and  plates  to  be  formulated  as  integral  equations,  in  terms  of  stress  intensity  factors 
at  the  crack  tip,  without  the  limitations  imposed  on  accuracy  by  beam  theory  approximations.  The 
proposed  functions  are  applied  to  investigate  the  influence  of  the  orthotropy  of  the  material  on  the 
fracture  behavior  of  double  cantilever  beams  in  the  presence  of  large  scale  bridging. 

Keywords:  Weight  function,  beam  theory,  double  cantilever  beam,  orthotropy  rescaling,  bridging 
mechanisms. 

1.  Introduction 

Over  the  last  decade,  an  interesting  class  of  fracture  mechanics  problems  has  arisen  in  laminated  beams 
and  plates,  in  which  nonlinear  crack  bridging  or  cohesive  mechanisms  act  over  several  scales.  One 
archetypal  material  system  of  particular  engineering  interest  is  laminated  polymer  composite  plies 
reinforced  through  the  thickness  by  stitches,  pins,  or  other  means.  Another  generic  problem  that  is 
important  in  electronics  and  packaging  is  the  separation  of  elastic  materials  connected  by  layers  of 
viscous  fluids  or  visco-plastic  media  such  as  solder.  In  all  these  delamination  problems,  crack  wake 
mechanisms  may  include  near  tip  processes  and  distinct  processes  acting  over  the  far  crack  wake.  In 
the  case  of  polymer  composites  containing  through-thickness  reinforcement,  the  near  tip  process  could 

To  whom  correspondence  should  be  addressed. 

Current  address:  Civil  and  Environmental  Engineering  Department,  Northwestern  University,  2145  Sheridan  Rd., 
Evanston,  IL  60208,  USA.  Email:  r-massabo@northwestem.edu ;  Phone  /  fax:  1  847  467  4105  /  1  847  491  4011; 
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be  the  formation  of  a  craze  zone,  which  is  usually  confined  to  a  domain  ~  1  mm  or  less;  while  the  far- 
wake  process  could  be  the  action  of  the  through-thickness  reinforcement,  which  can  extend  over  many 
tens  of  mm.  Since  the  far-wake  process  covers  a  domain  that  is  often  much  larger  than  the  laminate 
half-thickness,  h,  the  fracture  problem  can  be  analysed  accurately  using  an  appropriate  order  of  beam 
theory.  But  if  the  crack  is  short  (of  the  same  order  as  or  less  than  h),  then  beam  theory  fails  and  a  more 
accurate  method  is  needed  to  calculate  either  the  crack  tip  energy  release  rate  or  the  crack  tip  stress 
intensity  factors.  Greater  accuracy  is  also  required  when  addressing  a  near-tip  process  such  as  a  craze 
zone,  where,  even  for  a  long  crack,  beam  theory  gives  too  limited  a  description  of  near-tip  fields  to 
account  accurately  for  the  effects  of  the  zone  on  the  fracture  parameters. 

The  finite  element  method  (FEM)  offers  a  popular  and  convenient  numerical  alternative  to  beam 
theory.  However,  in  the  presence  of  plane  stress  or  plane  strain  conditions,  FEM  involves  the  solution 
of  fully  two-dimensional  problems;  and  can  present  some  tiresome  meshing  challenges,  especially 
when  the  boundaries  of  the  bridging  or  cohesive  zones  move  during  crack  evolution.  Moreover,  where 
the  bridging  tractions  depend  on  the  crack  displacement,  which  is  commonly  the  case,  the  nature  of 
singularities  in  the  near-tip  crack  displacement  can  be  difficult  to  resolve  properly  via  FEM.  This  has 
led  in  some  instances  to  inaccurate  results  appearing  in  the  literature. 

J 

Thus  motivation  exists  for  reducing  the  problem  of  a  bridged  crack  and  cohesive  zone  to  a  system  of 
integral  equations  by  means  of  weight  functions.  Weight  functions,  even  when  expressed  as 
engineering  approximations,  can  be  assured  to  have  the  correct  asymptotic  forms  and  therefore  offer 
the  greatest  control  in  numerical  methods  over  singularities  either  at  the  crack  tip  or  at  discontinuities 
in  the  bridging  mechanism  in  the  crack  wake.  For  special  cases,  the  integral  equation  formulation  leads 
to  closed  form  solutions.  Furthermore,  the  integral  equation  formulation  reduces  plane  problems  to 
one-dimensional  problems,  with  distance  from  the  crack  tip  the  only  spatial  variable.  This  allows  very 
rapid  scanning  of  large  quantities  of  parametric  problems. 

For  mode  I  fracture  problems,  the  required  weight  function  expresses  the  mode  I  crack  tip  stress 
intensity  factor  generated  by  a  pair  of  unit  point  forces  acting  normally  to  the  fracture  surfaces  at  any 
given  distance  from  the  crack  tip.  From  this  weight  function,  the  net  crack  tip  stress  intensity  factor 
due  to  the  combined  actions  of  an  applied  load  and  all  bridging  or  cohesive  mechanisms  over  all  scales 
can  be  computed  by  simple  superposition  and  the  crack  opening  profile,  by,  for  example,  an  application 
of  Castigliano’s  theorem  (e.g.,  [1,2]). 

For  special  problems  of  cracks  in  infinite  or  semi-infinite  bodies  or  cracks  in  long  strips  that  are 
propagating  across  the  strip,  the  definition  of  the  mode  I  weight  function  is  greatly  simplified  by  the 
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observation,  based  on  analytical  results  and  fairly  extensive  numerical  studies,  that  the  weight  function 
is  independent  or  almost  independent  of  die  degree  of  anisotropy  in  an  orthotropic  body  (provided  the 
crack  plane  is  aligned  with  a  plane  of  elastic  symmetry)  ([3,4].  Thus  weight  functions  evaluated  for 
isotropic  bodies  can  be  used  with  good  accuracy  in  more  general  cases.  The  only  appearance  of  the 
degree  of  anisotropy  is  in  the  elastic  constant  that  relates  the  crack  tip  stress  intensity  factor  to  the 
energy  release  rate. 

In  contrast,  in  plates  or  beams  such  as  a  double  cantilever  beam  (DCB)  specimen,  a  strong  dependence 
on  anisotropy  does  exist  and,  one  suspects  for  this  reason,  no  reasonably  simple  expression  for  the 
weight  function  for  an  orthotropic  DCB  has  appeared  in  the  literature.  Since  the  degree  of  anisotropy 
in  a  common  laminate  material  can  be  high,  especially  for  a  polymer  composite  but  also  for  structural 
ceramic  composites  or  layered  systems  in  electronics,  treating  the  effect  of  anisotropy  on  the  weight 
function  is  essential  to  establishing  an  accurate  integral  equation  formulation.  This  problem  is  taken  up 
here.  The  required  approximation  to  the  weight  function,  which  is  canonical  in  the  sense  that  it  retains 
all  the  asymptotic  features  and  the  essential  mathematical  form  of  the  exact  solution,  is  derived  from 
heuristic  consideration  of  the  weight  function  for  an  isotropic  DCB  specimen  and  the  implications  of 
orthotropy  rescaling.  An  approximate  weight  function  is  also  proposed  and  validated  numerically  for 
the  special  case  of  short  cracks  (much  less  than  the  laminate  thickness)  in  isotropic  and  degenerate 
orthotropic  DCB  specimens.  The  last  weight  function  will  permit  a  single  model  to  treat  the 
progression  from  the  initiation  of  a  bridged  crack  to  its  propagation  as  a  long  crack  with  large  scale 
bridging. 


2.  Mode  I  crack  tip  stress  intensity  factor  in  isotropic  double  cantilever  beams 

An  exact  solution  for  the  stress  intensity  factor  K&  due  to  a  pair  of  concentrated  forces  P  applied  per 
unit  width  onto  the  crack  faces  of  a  double  cantilever  beam  at  a  distance  d  from  the  crack  tip  has  been 
obtained  by  Foote  and  Buchwald  (see  Fig.  l.a)  [5].  They  solved  the  plane  elasticity  problem  by 
applying  the  Wiener-Hopf  technique  to  an  isotropic  infinite  strip  arbitrarily  loaded  by  opening  tractions 
acting  along  the  crack  surfaces  and  representing  the  concentrated  loads  in  terms  of  the  Dirac  delta 
function.  The  strip  and  the  applied  loads  were  assumed  to  be  symmetric  about  the  mid-plane.  A  simple 
formula  approximating  the  exact  solution,  which  has  an  accuracy  of  1.1%  and  can  be  applied  to  double 
cantilever  beams  with  an  uncracked  ligament  c  >  2h  and  containing  a  sufficiently  long  crack  is  also 
given  in  [5]: 
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where  h  is  the  half  thickness  of  the  specimen.  For  very  small  d/h  the  dimensionless  K\p  of  Eq.  (1) 
approaches  the  exact  solution  of  Irwin  for  a  semi-infinite  crack  in  an  infinite  sheet  [6]: 


For  very  large  d/h,  Eq.  (1)  approaches  the  elementary  beam  theory  solution  for  an  Euler-Bemoulli 
double  cantilever  beam  with  built-in  ends,  K\h°'s/P  =  -Jl2  d/h  [6].  For  d/h  >  0.3  the  exact  A^p  is  also 
well  represented,  with  an  error  always  lower  than  4%,  by  Gross  and  Srawley’s  boundary  collocation 
solution,  Kih°'5/P  =  -J\2  (d/h+ 0.687),  which  is  obviously  numerically  close  to  the  first  bracketed  term 
on  the  right  hand  side  of  Eq.  (1)  [6].  The  same  limit  solution  is  given  by  the  “augmented  beam  model” 
of  Kanninen  [7],  which  accounts  for  the  elasticity  of  the  uncracked  ligament  ahead  of  the  crack  tip  and 
removes  the  assumption  of  built-in  ends.  The  dimensionless  stress  intensity  factor  of  Eq.  (1)  is  shown 
in  the  semi-logarithmic  diagram  of  Fig.  l.b  along  with  the  two  limiting  solutions  for  small  and  large 
d/h. 

The  study  in  [5]  and  Eq.  (1)  refer  only  to  an  infinite  strip  or  a  double  cantilever  beam  with  an  infinitely 
long  crack  (the  crack  length  a  does  not  appears  in  Eq.  (1)).  A  lower  limit  for  the  normalized  crack 
length  a/h  of  the  double  cantilever  beam  specimen  must  be  set  for  Eq.  (1)  to  be  valid  for  all  d/h,  with  0 
<  d/h  <  a/h.  The  asymptotic  limit  of  Eq.  (1)  for  very  small  d/h,  given  by  Eq.  (2),  is  in  fact  correct  only 
if  a  »  d,  and  should  be  replaced  by  Tada’s  solution  [6]  for  a  finite  crack  of  length  a  in  a  semi-infinite 
sheet  when  a/h  also  becomes  very  small: 

Kjh °'5  =  2  Jh  1 
P  Jnci  ^l+x,/a 

where  x\  =  a  -  d  is  the  distance  of  the  applied  forces  from  the  crack  mouth.  The  value  of  K\  given  in 
Eq.  (3)  ranges  between  the  value  given  in  Eq.  (2),  for  a  »  d,  and  1.84  times  the  value  of  Eq.  (2),  for  a 
=  d. 

A  conservative  lower  limit  for  a/h  for  the  applicability  of  Eq.  (1)  can  be  set  as  a/h  =  0.3,  so  that  when 
d/h  =  a/h  -  0.3,  Gross  and  Srawley’s  solution  is  already  approached  by  the  exact  solution  and  the  error 
is  lower  than  4%.  This  lower  limit  has  been  validated  through  finite  element  calculations  which  also 
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1.3-0.3P-J  ,  (3) 
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showed  that  the  relative  error  between  the  predictions  of  Eq.  (1)  and  the  finite  element  results  is  within 
1.5%  for  all  normalized  crack  lengths  a/h  >  0.35  (see  Appendix  B). 


3.  Mode  I  crack  tip  stress  intensity  factor  in  orthotropic  double  cantilever  beams 

A  double  cantilever  beam  with  a  crack  of  length  a  loaded  by  a  pair  of  opening  forces  P\  acting  per  unit 
width  at  the  coordinate  xu  is  shown  in  Fig.  2.a.  Plane  stress  conditions  are  assumed  in  the  plane  X1-X3. 
The  material  is  linearly  elastic  and  orthotropic  with  principal  material  axes  xi  and  *3. 

As  shown  in  Appendix  A,  the  response  of  the  beam,  and  in  particular  the  stress  intensity  factor  at  the 
crack  tip  Kn,  is  controlled  by  two  dimensionless  material  parameters,  X  and  p,  which  depend  on  the 
four  elastic  constants  of  plane  orthotropic  elasticity.  In  terms  of  the  engineering  constants,  X  and  p  take 
the  form: 


P  = 


2G,3 


S' 


(4.a) 


where  E\  and  £3  are  the  Young’s  moduli  in  the  x\  and  X3  directions,  Gu  is  the  shear  modulus  and 
vl3  and  v3l  are  Poisson’s  ratios  in  the  plane  X1-X3.  The  parameter  n,  which  depends  on  p  and  will  be 

used  in  the  following,  is  also  introduced  in  (4.a).  To  ensure  the  positive  definiteness  of  the  strain 
energy  density  function,  the  orthotropy  ratios  are  such  that  X>  0  and  p  >  - 1.  In  typical  woods  and 
composites  p  is  in  the  range  0  <  p  <  5  and  X  in  the  range  0  <  X  <  20,  with  X  <  1  for  fiber  reinforced 
laminated  composites. 


For  plane  strain  conditions,  die  parameters  X  and  p  are  again  given  by  (4.a)  but  substituting  the  plane 
strain  constants: 


vlps 


A 


0-~Vl2V2 1) 


^3ps  ~ 


(I-VjjVjj) 


’■'Bps  “ 


(l-vnv21)  ’ 


(^31  +  ^21)  ? 

(I-V32V23) 


(4.b) 


for  Ei,  £3 ,  v13  and  v31 . 

In  the  case  of  isotropy,  X  =  p  =  1  and  the  stress  intensity  factor  at  the  crack  tip  of  the  specimen  of  Fig. 
2  is  given  by  Eq.  (1),  with  the  substitutions  d-a-x u  and  P  =  P\.  The  stress  intensity  factor  for  the  case 
of  orthotropy  will  be  derived  in  the  following. 
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3.1  Solution  for  orthotropic  double  cantilever  beams  with  large  notches 

An  expression  for  the  strain  energy  release  rate  $  in  the  orthotropic  double  cantilever  beam  specimen 
of  Fig.  2.a  has  been  obtained  by  Suo  et  al.  [8]  making  use  of  the  orthotropy  rescaling  technique  and 
numerical  calculations.  The  expression  of  $  given  in  [8]  depends  on  two  terms:  the  first  is  the 
elementary  beam  theory  result,  which  does  not  depend  on  the  orthotropy  and  is  an  exact  elasticity 
asymptote  as  (a-x\\)/h  oo;  the  second  is  a  correction  term  which  accounts  for  the  influence  of  finite 
(a-x\[)/h.  The  expression  proposed  in  [8]  has  an  accuracy  of  1%  for  all  (a-xn)/h  £  2Xm  and  0<p<4. 

The  stress  intensity  factor  at  the  crack  tip  can  be  obtained  from  £i  of  [8]  using  the  orthotropic 
relationship  between  and  K\\ 


where 


•^i  —  > 


(5.a) 


£,A3'4 

V(l  +  P)/2 


£,/l3/4 

n 


(5.b) 


In  a  dimensionless  form,  Ku  is  then  given  by: 


KKh01 

P> 


h 


1+Y, \{p)X 


-1/4 


where: 

Y \{p)  =  0.677  +  0.146(p  - 1)  -  0.0178(p  - 1)2  +  0.00242(p  - 1)3  ,  (7) 

is  a  correction  term  that  accounts  for  the  influence  of  p  and  has  been  obtained  in  [8]  through  finite 
element  calculations.  For  an  isotropic  material,  X  =  p  =  1,  Yi(p)  =  0.677  and  Eq.  (6)  coincides  with 
Gross  and  Srawley’s  solution.  The  last  term  on  the  right  hand  side  of  Eq.  (6)  would  not  be  present  if  the 
two  delaminated  beams  were  assumed  to  be  Euler-Bemoulli  beams  built  in  at  the  crack  tip.  It  describes 
the  influence  of  the  elasticity  of  the  uncracked  ligament  ahead  of  the  crack  tip.  It  vanishes  when  the 
ratio  {a-x\\)/h  becomes  very  large,  because  in  this  limit  the  condition  of  built-in  ends  becomes 
increasingly  well  satisfied. 
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Equations  (5)  and  (6)  correctly  define  the  stress  intensity  factor  and  the  strain  energy  release  rate  in 
conventional  2D  laminates  tested  with  typical  notched  specimens,  for  which  (a-xi\)/h  >  2A'U4.  On  the 
other  hand,  Eq.  (6)  cannot  be  applied  to  investigate  the  influence  of  opening  or  closing  forces  acting 
very  close  to  the  crack  tip  such  as  those  produced  by  the  bridging  mechanisms  due  to  through-thickness 
reinforcement  or  crazing.  To  analyze  these  cases  a  solution  for  J^i  valid  over  the  whole  range  {a-x\\)/h  > 
0  is  required. 

3.2  Solution  for  degenerate  orthotropic  double  cantilever  beams 

Consider  the  rescaled  problem  of  Fig.  2.b  where  the  xi  axis  of  the  problem  of  Fig.  2.a  is  rescaled  as  £= 
A1'4  x\.  Following  the  orthotropy  rescaling  technique  recalled  in  the  Appendix  A,  the  geometry  and 
boundary  conditions  are  rescaled  in  such  a  way  that  the  crack  length  becomes  Ama,  the  forces  applied 
per  unit  width  XV4P\  and  the  stress  intensity  factor  at  the  crack  tip  XViK\\. 


As  shown  in  Appendix  A,  in  a  degenerate  orthotropic  solid  characterized  by  p  =  1,  foe  rescaled 
problem  in  the  plane  £-*3  is  governed  by  foe  same  equations  of  foe  original  problem  in  an  isotropic 
solid.  Consequently,  from  Eq.  (1),  the  rescaled  stress  intensity  factor  takes  foe  form: 


where  the  term  0.673  of  Eq.  (1)  has  been  changed  to  0.677  to  satisfy  Eqs.  (6)  and  (7)  for  large  {a-x\\)/h. 
In  a  dimensionless  form,  similar  to  that  of  Eq.  (1),  Ky,  becomes: 


Equation  (9)  is  an  exact  extension  of  foe  isotropic  result  (1)  to  a  degenerate  orfootropic  material.  The 
equation  has  been  validated  through  finite  element  calculations  for  A  varying  in  foe  range  0.025  <,  A  <, 
1.0  and  long  enough  cracks  with  error  always  within  2%  (see  Appendix  B).  The  equation  applies  to 
beams  with  an  uncracked  ligament  c/h  >  2XU4. 
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Equation  (9)  includes  the  result  expected  for  very  small  (a-x\\}/h,  tending  asymptotically  to  Irwin’s 
solution  of  Eq.  (2).  In  this  limit,  the  stress  intensity  factor  is  not  affected  by  the  anisotropy  of  the 
material,  as  demonstrated  by  Sih  et  al.  [3]  for  cracks  in  infinite  sheets  loaded  by  self-equilibrating  loads 
acting  on  the  crack  surfaces.  For  large  (a-x\\)/h,  Eq.  (9)  approaches  the  solution  of  Suo  et  al.,  Eq.  (6), 
with  p  =  n  =  1  and  Y\{p)  =  0.677.  Note  that  this  upper  limit  as  well  as  the  connecting  function  given 
by  the  third  term  on  the  right  hand  side  of  Eq.  (9)  depend  on  the  orthotropy  of  the  material  through 
A.  The  dimensionless  stress  intensity  factor  of  Eq.  (9)  is  depicted  in  Fig.  3  for  different  values  of 
X.  The  dotted  curves  in  the  figure  show  the  limiting  solutions. 

Equation  (9)  is  correct  for  all  {a-x\\)/h  provided  the  crack  length  a/h  is  higher  than  a  limit  value  which 
the  rules  of  orthotropy  rescaling  sets  equal  to  a/h  =  0.35A"I/4  for  a  1.5%  accuracy.  Finite  element 
calculations  show  that  this  is  in  fact  a  conservative  limit  for  all  cases  with  A  <  1.0  (see  Appendix  B). 


3.3  Solution  for  orthotropic  double  cantilever  beams 


The  stress  intensity  factor  in  generally  orthotropic  (A  ^  1  and  p  &  1)  double  cantilever  beams  must 
approach  the  known  asymptotic  solutions  for  large  and  small  {a-x\\)/h.  For  very  small  ratios  (a-x\\)/h 
the  limit  solution  is  given  by  Irwin’s  result  for  a  semi-infinite  crack  in  an  infinite  sheet,  Eq.  (2).  As 
already  noted,  K\  in  this  limit  does  not  depend  on  the  elastic  constants.  For  large  {a-x\\)/h,  the  limit 
solution  is  given  by  Eq.  (6).  Solutions  for  intermediate  values  of  {a-x\\)/h  require  rigorous  analyses  of 
the  problem  which  have  been  performed  through  finite  element  calculations.  Only  the  effect  of  p  needs 
to  be  calibrated  numerically  since  the  dependence  of  A  is  known  analytically  and  in  materials  with 
cubic  symmetry,  A  =  1,  the  plane  elasticity  problem  is  controlled  by  p  alone  (see  Appendix  A). 

Based  on  these  observations,  the  following  expression  is  proposed  as  a  canonical  approximation  to  Kn 
for  the  degenerate  orthotropic  material,  derived  by  modifying  Eq.  (9)  to  have  the  right  asymptotic 
behavior  and  fit  the  finite  element  results: 


0.5  ^3/8 
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(10) 


The  validity  of  Eq.  (10)  has  been  checked  through  finite  element  calculations  for  X  and  p  in  the  range 
0.025  <X<  l  and  1  < p <  5  with  errors  always  less  than  2%.  Equation  (9)  is  applicable  for  all  ratios  (a- 
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x\\)/h  provided  the  crack  length  a/h  is  higher  than  a  limit  value;  the  limit  previously  defined  for  the 
degenerate  orthotropic  materials,  a/h  =  0.35  A'1/4,  proves  to  be  a  conservative  limit  for  all  values  of 
p( see  Appendix  B).  As  for  the  degenerate  orthotropic  material,  Eq.  (10)  applies  only  to  beams  with 
uncracked  ligaments  cfh  >  2A',/4. 

The  dimensionless  stress  intensity  factor  of  Eq.  (10)  is  shown  in  Fig.  4  for  A  =  1  and  different  values  of 
the  orthotropy  parameter  p  to  highlight  the  influence  of  p.  For  large  ratios  (a-x\\)/h,  the  effect  of 
p  becomes  increasingly  strong,  as  expected  from  Eq.  (6),  and  the  relative  error  between  the  exact 
solution  and  the  solution  for  the  degenerate  orthotropic  material  becomes  (1-1/Vw ).  Thus  only  for 
rather  small  cracks,  with  (a-xu)/h  <  0.7,  could  the  solution  of  Eq.  (9)  for  the  degenerate  orthotropic 
solid  (p=  1)  be  used  for  all  values  of  p  with  acceptable  error.  For  general  cracks,  Eq.  (10)  must  be 
used.  The  dependence  of  the  weight  function  on  p  for  larger  cracks  is  not  shared  by  the  strain  energy 
release  rate,  £i,  which  becomes  independent  of  p  and  A  when  (a-x\\)/h  is  large  (see  Eqs.  (5-6)). 

Figure  5  shows  the  dimensionless  stress  intensity  factor  of  Eq.  (10)  for  three  pairs  of  A  and 
p  representative  of  composites  used  in  the  aeronautical  industry,  A=  0.1  and  p  =  3  (e.g.  a  quasi¬ 
isotropic  carbon-epoxy  laminate),  A  =  0.05  and  p  =  5  (e.g.  a  unidirectional  boron-epoxy  laminate)  and 
A  =  1  and  p  -  1  (an  isotropic  material). 

4.  Mode  I  crack  tip  stress  intensity  factor  in  double  cantilever  beams  with  short  cracks 

The  stress  intensity  factors  obtained  in  the  previous  sections  refer  to  DCB  specimens  with  normalized 
lengths  of  the  crack  higher  than  a/h  =  0.35A'1/4.  To  the  authors’  knowledge  there  are  no  expressions 
available  in  the  literature  for  cracks  of  lengths  a/h  <  0.35A‘1/4,  either  for  isotropic  or  orthotropic 
specimens.  Knowing  the  expression  of  the  stress  intensity  factors  for  cracks  of  any  length  can  be  useful 
in  dealing  with  various  mechanical  problems,  such  as  the  problem  of  crack  initiation. 

In  the  case  of  isotropy,  the  problem  is  now  controlled  by  two  dimensionless  variables,  which  define  the 
normalized  distance  of  the  applied  forces  from  the  crack  tip,  (a-x\\)/h,  and  their  relative  position  with 
respect  to  die  crack  length  x\\/ a.  This  behavior  can  be  observed  in  Eq.  (3),  which  gives  the  asymptotic 
limit  of  the  stress  intensity  factor  for  very  small  (a-x\\)/h.  Equation  (1)  gives  the  exact  stress  intensity 
factor  for  very  long  cracks  and  it  is  therefore  the  asymptotic  limit  for  x\  \/  a  -»  1  of  the  solution  being 
sought.  Equation  (1)  can  then  be  modified  as  follows: 
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where  the  second  term  on  the  right  hand  side  has  been  changed  from  that  given  by  Eq.  (2)  to  that  given 
by  Eq.  (3)  and  the  multiplying  function  fix\\/a,  {a-x\\)/h)  has  been  introduced  in  the  third  term.  The 
functional i/o,  {a-x\\)/h)  should  be  unity  for  x\\/  a  ->  1,  give  the  right  asymptotic  behaviors  for  large 
and  small  ratios  {a-x\\)/h,  and  concur  with  numerical  solutions  for  intermediate  values  of  {a-x\\)/h  and 
xu/a  <  1.  The  following  is  such  a  function: 


(12.a) 


(12.b) 


(12.c) 


The  particular  shape  of  the  functional  i/o,  (o-X[\)/h)  has  been  chosen  so  that/is  almost  equivalent  to  g 
over  most  of  the  range  0  <xu/a  <  1;  and  only  forxn/ a  — »  1  becomes  equal  to  1.  Equation  (11)  has  an 
accuracy  better  than  2%  for  xy,/ a  <  1  and  (a-x\  \)/h  >  0.  It  applies  to  beams  with  an  uncracked  ligament 
c/h>  2. 


Figure  6  shows  the  dimensionless  stress  intensity  factor  of  Eq.  (1 1)  as  a  function  of  the  normalized 
distance  of  the  applied  forces  from  the  crack  tip,  {a-x\-yh,  on  varying  x\\/  a.  The  figure  highlights  the 
expected  variations  in  Ku  with  respect  to  the  solution  of  Eq.  (1),  shown  in  the  diagram  as  the  curve 
marked  xy,/a=  1. 
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Following  the  procedure  of  Section  3.2,  Eqs.  (11)-(12)  can  be  rescaled  to  define  the  mode  I  crack  tip 
stress  intensity  factor  in  degenerate  orthotropic  DCB  specimens: 


Equation  (13)  has  been  validated  through  finite  element  calculations  for  X  varying  in  the  range  0.025  < 
A<  1.0  with  error  always  less  than  2%.  It  applies  to  beams  with  an  uncracked  ligament  c/h  >  2X'm.  The 
influence  of  X  on  the  solution  is  similar  to  that  already  observed  in  Figure  3.  The  extension  of  Eq.  (13) 
to  a  generally  orthotropic  material  is  more  complicated. 

5.  Summary  of  mode  I  weight  functions  for  orthotropic  double  cantilever  beams 

The  stress  intensity  factors  calculated  in  the  previous  section  define  the  mode  I  weight  functions  for 
double  cantilever  beams  with  a  mid-plane  crack  of  length  a,  hi(xu,a,h,p,X)  =  K\\  l  P\,  where  Ku  is  given 
by  Eq.  (1)  and  Eq.  (11)  for  isotropic  materials,  by  Eq.  (9)  and  Eq.  (13)  for  degenerate  orthotropic  and 
by  Eq.  (10)  for  generally  orthotropic  materials.  The  results  are  summarized  in  a  dimensionless  form  in 
Table  1  and  Table  2. 

The  weight  functions  allow  the  definition  of  the  stress  intensity  factor,  Kip,  at  the  tip  of  the  crack  in 
double  cantilever  beams  subject  to  general  distributions  of  tractions,  p(xi),  acting  along  the  crack  faces: 


ll 


(15) 


KiP=  jfJi(xl,a,h>pJ)p(xl)dxl  , 
a0 

where  ao  and  xif  are  the  limits  of  the  loaded  region.  In  large  scale  bridging  problems  ao  typically 
represents  the  length  of  a  notch  or  a  region  of  the  laminate  where  there  are  no  active  bridging 
mechanisms  and  xk  is  the  tip  of  the  crack  (Fig.  7). 


6.  Fracture  behavior  of  orthotropic  double  cantilever  beams  in  the  presence  of  large  scale 
bridging 

The  weight  functions  derived  above  are  applied  in  this  Section  to  analyze  crack  growth  in  an 
orthotropic  double  cantilever  beam  with  bridging  mechanisms  acting  along  the  crack.  The  stress 
intensity  factor  at  the  crack  tip  of  the  beam  when  loaded  by  a  pair  of  opening  forces  P,  acting  at  the 
coordinate  *i  =  0,  and  by  closing  tractions  p,  acting  along  a-ao  as  shown  in  Fig.  7,  is  given  by: 

afh 

Ki=Kir+Klf  =  Klf-  / h,a/ h,p,Z)h05p(xl )  d(-*-)  ,  (16) 

C^th  " 

where  Ajp  is  given  by  Eq.  (10)  for  xn  =  0,  and  K\p  is  the  stress  intensity  factor  due  to  opening  tractions 

P- 

If  the  tractions  p  are  created  by  through-thickness  reinforcement  or  other  bridging  mechanisms  acting 
along  the  crack  surfaces,  they  depend  on  the  crack  opening  displacement,  w3,  and  are  a  priori  unknown. 
In  particular,  if  w3(x  1)  >  0,  then  p[u3(xx )]  =  p3[u3(x  1)]  are  closing  tractions  depicting  the  bridging 
mechanisms.  The  value  of  p3  as  a  function  of  w3  is  defined  through  the  bridging  traction  law,  p3{u3\ 
which  is  deducible  from  micromechanical  or  macromechanical  models  and  experiments  (see  [9]  for  the 
case  of  stitched  laminates).  If  u3{x\)  =  0,  then  p[u3(x\)]  =  -  pc[(*i)]  are  opening  tractions  depicting  the 
effect  of  contact  pressure  arising  along  the  crack.1  The  contact  pressure  and  the  size  of  the  regions  of 
contact  or  bridging  are  unknown  a  priori  and  can  be  determined  through  a  compatibility  condition  for 
the  crack  opening  displacement. 

At  the  onset  of  crack  propagation,  K\  of  Eq.  (17)  is  equal  to  the  intrinsic  fracture  toughness,  K\  =  K\c, 
and  the  dimensionless  critical  load  for  crack  propagation  takes  the  form: 


1  The  existence  of  regions  of  contact  along  the  surfaces  of  the  crack  in  double  cantilever  beams  in  large  scale  bridging 
conditions  has  been  observed  experimentally  in  [10,1 1]  and  explained  theoretically  in  [12]. 
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where  pw  is  a  normalizing  value  of  the  crack  face  tractions,  pi,  given  for  instance  by  their  maximum 
value.  Recalling  the  expression  for  E\  and  that  K\c  =  ,  Eq.  (18)  can  be  modified  to  facilitate 

direct  comparison  between  isotropic  and  orthotropic  cases: 


P*  _=  1  U  I  Pso!i_  °f  h(x  /ha/h  a  juOiPWVldfi)  (\g\ 

J&kh  b  (0,  a/h,p,  X)h0i  |  V»  J1‘  ‘  ’  P,o 


As  we  will  see  later,  the  dimensionless  number  appearing  within  the  right  hand  side  of  Eq.  (18)  and 
(19),  p30A0'5  or  p30h05A31*  !{Klc4n) ,  is  a  measure  of  the  brittleness  of  the  structure  [13]. 


The  crack  opening  displacement  at  the  generic  coordinate  xn  is  defined  as  a  function  of  the  crack  tip 
stress  intensity  factor  and  the  elastic  constants  of  the  material  through  the  following  relationship: 

a 

»3(»ii)  =  «j(*n>p+'«s(*ii)p  “  V*-  (20) 

«0 

where  Ajj  is  the  localized  compliance,  which  defines  the  crack  opening  displacement  at  xn  due  to  a  pair 
of  unit  opening  forces  acting  at  xij.  The  localized  compliance  can  be  defined  through  an  energy  balance 
or  Castigliano’s  theorem  and  is  given  by: 

=  17 )b (*u  / b a'/ h, p, A)  b (x,j  / h,a't h, p, A) da'  (21) 

E\  # 

where  Ex  is  the  orthotropy  constant  of  Eq.  (5.b). 


Once  the  crack  opening  displacement  is  known,  the  statically  indeterminate  problem  defined  by  the 
nonlinear  integral  equations  (19),  (20)  and  (21)  is  solved  for  general  bridging  laws,  ps{ui),  through  a 
discretization.  A  self-consistent  solution  for  the  crack  profile  is  obtained  iteratively  through  a 
numerical  procedure  [1,2]. 

In  the  special  case  of  bridging  mechanisms  described  by  a  Dugdale  type  bridging  law,  pz  -  pzo,  the 
bridging  tractions  do  not  depend  on  u$  and  Eq.  (19)  alone  gives  the  dimensionless  critical  load  for 
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cmck  propagation.  The  model  predicts  fite  absence  of  regions  of  contact  for  this  case  confirming  the 
beam  theory  predictions  of  [  1 2] . 

Figures  8  a  8  b  and  8.c  show  dimensionless  diagrems  of  the  critical  load  for  crack  propagation  as  a 
of  the  dimensionless  parameter  pi0h  ic  1  ,  , 

P  3  ^bteined  k°[12]  assuming  the  two  delaminated  beams  to  be  built  in  at  the  crack  tip  and 
—  s  U  -  -so  shown  in  fite  diagrams.  The  infiuence  of  shear 

deformations  (y*  0)  on  the  beam  theory  results  is  shown  in  Figure  8.b. 

The  response  of  fite  double  cantilever  beam  is t  — 1^  ^^t  S 
n  h^uEi:  and  the  two  orthotropy  ratios,  A  and  p.  The  three  aiagram  * 

trisition  in'flte  macrestrecmra,  response  of  fite  beam  which  is  typical  »f 

Dugdale-^e  reinforcement  113]:  fite  fa  are  changed  so  as 

increasing  the  dimensionless  group  piah  /^i  •  Moreover,  i  ^30,  yvu 

■  1  r'kanctnnt  the  beams  will  experience  the  same  type  of  response.  As  for 

8.a)  and  facilitates  strain  hardening  behavior  (Fig.  8.b-c). 

T*  anisotropy  of  fite  materia,  affects  fite  response  —  ""'a 

case  of  largpcale^n^"g  corresponding  to  the  point  where  all  curves  cross 

transition  value  of  a/h  1/(1.  /ip30“  yyic  i  , 

each  ofiter.  ,f  *  is  smaller  than  fite  fiansirion  value,  even 

common  slope  given  by  pi0h  /( )»  ana  me  u  1/ov.  ...m 

solufion  andthe  correct  solution  b^njj-  ~td" of^taLic  Lire 
P^!&-  T*  fractional  ereor  »  Lons  and  coincides  wifi,  the 

toughness  of  the  laminate,  an  ^  181  It  depends  only  on  the  crack  length  and  the 

Z  loved  finough  fite  use  of  modified  beam  decries 

accounting  for  the  influence  of  file  elastic  material  ahead  of  the  crack  t,p  [7]. 
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7.  Conclusions 


Canonical  expressions  for  the  mode  I  weight  function  for  an  orthotropic  DCB  specimen  have  been 
derived  and  validated  by  finite  element  methods.  By  combining  results  for  short  and  long  cracks,  a  set 
of  expressions  has  been  assembled  that  are  collectively  accurate  to  within  2%  over  a  wide  range  of 
degrees  of  anisotropy.  The  canonical  expressions  retain  the  correct  asymptotic  limits  for  small  and 
large  cracks  and  special  cases  of  elastic  symmetry.  Thus  expedient  and  accurate  study  of  singularities 
and  other  special  characteristics  of  bridged  crack  problems  can  be  based  on  their  use. 

The  weight  functions  are  the  essential  prerequisite  for  integral  equation  formulations  of  fracture 
problems  in  which  bridging  or  cohesive  mechanisms  are  present  on  small  and  large  scales.  The  weight 
functions  will  allow  equally  accurate  and  convenient  solution  of  crack  initiation  problems,  in  which  the 
crack  length  remains  smaller  than  or  comparable  to  the  laminate  half-width;  and  large-scale  bridging 
problems,  in  which  bridging  effects  may  extend  over  zones  many  tens  of  times  the  laminate  half-width. 
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Appendix  A  -  Orthotropy  rescaling 


An  orthotropy  rescaling  technique  has  been  proposed  by  Suo  [14]  to  reduce  plane  elasticity  problems 
for  orthotropic  materials  to  equivalent  problems  for  materials  with  cubic  symmetry. 


For  plane  stress  problems  in  the  plane  X1-X3,  with  jci  and  *3  principal  material  axes  and  assuming 
negligible  body  forces,  compatibility  implies: 


etc,4  dx^dx/ 


=  0 


(A.1) 


The  dimensionless  parameters  A  and  p  in  Eq.  (A.l)  depend  on  the  elastic  constants  of  the  material  and 
have  been  defined  in  the  body  of  the  text,  Eq.  (4.a). 


Equation  (A.l)  shows  that  in  orthotropic  simply-connected  sheets  with  tractions  prescribed  on  the 
boundary,  stresses  depend  on  only  two  materials  parameters,  A  and  p.  Moreover,  if  the  xi  axis  is 
rescaled  as  £  =  Am  x\,  Eq.  (A.  1)  shows  dependency  on  p  only: 


d4U  _  d4U  d4U 
— r+2/> — n — ;r+ — r 
d%4  d?dx32  dx34 


(A.3) 


Equations  (A.l),  (A.3)  and  (A.4)  show  that  the  orthotropic  problem  in  the  rescaled  plane  £-*3  is 
reduced  to  the  equivalent  problem  for  a  material  with  cubic  symmetry,  namely  with  A  -  1.  Moreover, 
the  special  case  of  a  degenerate  orthotropic  material,  a  material  with  p  =  1,  is  governed  in  the  rescaled 
plane  by  the  same  equation  which  governs  the  isotropic  problem.  Consequently,  solutions  for 
degenerate  orthotropic  materials  can  be  readily  constructed  from  the  isotropic  results. 
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Equations  (A.3)  and  (A.4)  also  show  that  the  stress  intensity  factors  for  a  crack  in  the  *i  direction,  ft 
and  Km  being  given  by  the  limits  for  x,->0  of  (2kx,)’%(*i>0)  and  (2nxl)l/2a13(x1,0),  respectively, 
rescale  as  Xv%  and  Xmft\.  Moreover,  the  component  of  the  resultant  force  on  an  arc  2*  rescale  as  X 
v%  while  the  component  T\  remains  unchanged. 
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Appendix  B  -  Finite  element  validation 


Finite  element  calculations  were  performed  using  ANSYS5.5  [15]  in  order  to  validate  the  proposed 
weight  functions  and  their  limit  of  applicability.  The  stress  intensity  factor  generated  by  a  pair  of  point 
forces  acting  along  the  crack  surfaces  has  been  calculated  using  two  different  techniques,  the 
displacement  correlation  technique  and  the  J-integral.  The  finite  element  meshes  were  composed  of 
quadratic  elements  with  a  rosette  of  quarter-point  elements  at  the  crack  tip  that  embody  the  linear 
elastic  fracture  mechanics  square  root  singularity.  Convergence  studies  were  carried  out  for  each 
loading  configuration  to  define  the  necessary  degree  of  refinement  of  the  meshes. 


The  displacement  correlation  technique  used  here  has  been  proposed  by  Shih  et  al.  [16]  on  the  basis  of 
Barsoum’s  demonstration  that  quarter-point  elements  show  the  square  root  singularity  of  elastic 
fracture  mechanics  [17].  The  technique  has  then  been  extended  to  anisotropic  problems  by  Sauma  and 
Sikiotis  [18].  It  allows  the  calculation  of  stress  intensity  factors  performing  simple  post-processor 
operations  on  4  nodal  displacements  of  the  quarter-point  elements  along  the  crack  faces.  The  J-integral 
calculations  [19]  were  performed  following  the  procedure  implemented  in  ANSYS5.5.  The  stress 
intensity  factor  was  then  extracted  from  7=  $  using  the  relationship  (5.a).  Due  to  the  very  fine  meshes 
used  in  the  analyses  the  two  calculations  yielded  virtually  identical  results. 

Validation  of  Eqs.  (9)  and  (10)  for  large  cracks 

Table  I  shows  normalized  stress  intensity  factors  and  relative  errors  between  FE  calculations  and  Eq. 
(9)  for  a  degenerate  orthotropic  material  (p=  1)  with  different  values  of  X  in  the  range  0.025  <  X  <  1.0 
on  varying  (a-X\{)/h  in  the  range  0.01  -  5.0. 


TABLE  I 


p=  1 

5.0 

2.0 

1.0 

0.5 

0.3 

0.1 

0.05  I 

0.01 

19.666 

4.112 

3.605 

4.314 

Bm 

19.643 

mmm 

4.161 

3.608 

4.300 

umm 

Error  (%) 

0.116 

0.342 

-1.318 

-0.085 

0.309 

0.682 

X  =  0.25 

FEM 

12.271 

6.092 

4.019 

3.103 

2.879 

4.090 

8.208 

Eq.(9) 

12.243 

6.070 

4.054 

3.146 

2.908 

0.348 

-0.871 

-1.384 

-1.009 

0.868 

bebm 

FEM 

5.822 

3.216 

2.382 

2.152 

2.214 

2.939 

3.859 

8.107 

5.802 

3.224 

2.434 

2.173 

2.218 

2.915 

3.825 

8.074 

Error  (%) 

0.343 

-0.247 

-2.151 

-0.977 

-0.179 

0.819 

0.858 

0.399 

Table  II  shows  normalized  stress  intensity  factors  and  relative  errors  between  FE  calculations  and  Eq. 
(10)  for  a  material  with  cubic  symmetry  (X=  1)  with  different  values  of  p  in  the  range  1  <  p<  5  on 
varying  {a-x\\)/h  in  the  range  0.01  -  5.0. 
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TABLE  II 


A=1 

(a-xvyh 

1.0 

0.5 

0.3 

0.1 

0.05 

0.01 

1 

See  Table  I 

n 

FEM 

17.237 

4.145 

3.635 

_ 

8.319 

Eq.(10) 

17.178 

8.408 

5.494 

4.099 

3.685 

4.353 

0.337 

1.057 

1.704 

1.118 

0.009 

-0.629 

-0.100 

1 

H 

FEM 

mm 

WESSM 

5.566 

4.247 

3.750 

3.736 

4.400 

8.341 

Eq.  (10) 

■ 

■jMjH 

3.730 

3.784 

4.435 

8.327 

Error  (%) 

0.555 

1.725 

2.172 

2.00 

0.517 

-1.261 

-0.788 

0.170 

Finally,  Table  III  shows  normalized  stress  intensity  factors  and  relative  errors  between  FE  calculations 
and  Eq.  (10)  for  a  generally  orthotropic  material  with  different  values  of  A  and  p  in  the  range  0.025  < 
A<  1.0  and  1  <  p  <  5  on  varying  {a-xx\)/h  in  the  range  0.01  -  5.0. 


TABLE  in 


(a-x^/h 

5. 0/6.0 

2.0 

1.0 

0.5 

0.3 

0.1 

0.05 

0.01 

1 

8.145 

4.460 

wmm 

2.608 

3.146 

4.001 

8.169 

Eq.  (10) 

8.086 

4.390 

WEEM 

2.627 

3.152 

3.991 

8.142 

Error  (%) 

0.734 

1.564 

0.71 

-0.186 

-0.724 

-0.206 

KE39I 

■tltlcV 

FEM 

6.327 

3.725 

2.841 

2.490 

2.472 

3.077 

KEB 

mmm 

6.235 

3.646 

WBm 

3.101 

3.961 

8.132 

1.447 

2.124 

-0.454 

-1.163 

-0.781 

-0.177 

0.190 

X  =  0.025 
p=5 

FEM 

5.165/5.837 

2.254 

2.299 

3.891 

8.118 

Eq.  (10) 

5.087/5.753 

3.095 

2.475 

2.275 

2.330 

3.005 

3.895 

8.106 

Error  (%) 

1.508/1.433 

1.202 

0.0892 

-0.942 

-1.306 

-0.705 

-0.095 

0.151 

Validation  of  Eqs.  (1),  (9)  and  (10)  for  short  cracks 

To  validate  Eqs.  (1),  (9)  and  (10)  for  very  short  cracks  and  define  the  minimum  crack  length  required 
for  their  applicability,  the  problem  of  a  double  cantilever  beam  loaded  by  a  pair  of  concentrated  forces 
acting  at  the  coordinate  jcu  =  0  has  been  considered  and  analyzed  by  the  fmite  element  method  on 
varying  the  orthotropy  ratios  in  the  ranges  0.025  <  A  <  1.0  and  1  <  p  <  5. 

For  an  isotropic  material,  Eq.  (1)  proves  to  have  an  accuracy  higher  than  1.5%  for  all  crack  lengths 
higher  than  a/h  =  0.35.  For  a  degenerate  orthotropic  material,  Eq.  (9)  proves  to  have  an  accuracy 
higher  than  2%  for  all  crack  lengths  higher  than  a/h  =  0.35A*1'4.  The  same  limit  is  also  a  conservative 
limit  for  generally  orthotropic  materials. 
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TABLE  CAPTIONS 


Table  1 :  Dimensionless  mode  I  weight  functions  for  isotropic,  degenerate  orthotropic  and  generally 
orthotropic  double  cantilever  beams.  Accuracy  better  than  2%  for  crack  lengths  a/h  > 
0.35A'1/4. 

Table  2:  Dimensionless  mode  I  weight  functions  for  isotropic  and  degenerate  orthotropic  double 
cantilever  beams.  Accuracy  better  than  2%  for  all  crack  lengths. 


FIGURE  CAPTIONS 

Figure  1  .a:  Schematic  of  an  infinitely  long  strip  [5] . 

Figure  l.b:  Solid  curve:  dimensionless  stress  intensity  factor  in  the  isotropic  strip  of  Fig.  l.a  (Foote  and 
Buchwald  and  Eq.  (1)).  Dashed  curves:  limiting  solutions  for  very  small  d/h  (Irwin)  and 
large  d/h  (Gross  and  Srawley). 

Figure  2:  a)  Schematic  of  a  double  cantilever  beam  loaded  by  a  pair  of  point  forces  P,.  b)  Orthotropy 
rescaling:  the  geometry  and  boundary  conditions  are  rescaled  from  the  original  problem  a). 

Figure  3:  Dimensionless  stress  intensity  factor  in  the  degenerate  orthotropic  (/? =  1)  beam  of  Fig.  2.a , 
Eq.(9). 

Figure  4:  Dimensionless  stress  intensity  factor  in  the  orthotropic  beam  of  Fig.  2.a  with  cubic 
symmetry  (X  =  1),  Eq.  (10). 

Figure  5:  Dimensionless  stress  intensity  factor  in  the  orthotropic  beam  of  Fig.  2.a,  Eq.  (10). 

Figure  6:  Dimensionless  stress  intensity  factor  in  the  isotropic  beam  of  Fig.  2.a  on  varying  xy,/ a,  Eq. 

(11). 

Figure  7:  Schematic  of  a  double  cantilever  beam  in  the  presence  of  bridging  mechanisms. 

Figure  8:  Dimensionless  critical  load  versus  normalized  crack  length  in  orthotropic  DCB  specimens: 

a)  No  bridging  mechanisms,  pi0h05  / ^ICE,  =  0.0;  b)  p30hos  / ■^QlcEl  =  0.1;  c) 

p30h0S 1 4$ ic^i  =  0-3.  (curves  (a):  X  =  1  and p-  1;  curves  (b)  X-  0.1  and  p=  3;  curves  (c) 

/l  =  0.05  and/?  =5). 
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Figure  l.b:  Solid  curve:  dimensionless  stress  intensity  factor  in  the  isotropic  strip  of  Fig.  l.a  (Foote  and 
Buchwald  and  Eq.  (1)).  Dashed  curves:  limiting  solutions  for  very  small  d/h  (Irwin)  and  large  d/h 

(Gross  and  Srawley). 
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Figure  2:  a)  Schematic  of  a  double  cantilever  beam  loaded  by  a  pair  of  opening  forces  P\. 
b)  Orthotropy  rescaling:  the  geometry  and  boundary  conditions  are  rescaled 
from  the  original  problem. 
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Figure  4:  Dimensionless  stress  intensity  factor  in  the  orthotropic 
beam  of  Fig.  2.a  with  cubic  symmetry  (A=  1),  Eq.  (10). 
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Figure  5:  Dimensionless  stress  intensity  factor  in  the  orthotropic  beam  of  Fig.  2.a,  Eq.  (10). 


29 


Figure  6:  Dimensionless  stress  intensity  factor  in  the  isotropic  beam  of  Fig.  2.a 
on  varying  xn  /  a,  Eq.  (11). 
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Figure  7:  Schematic  of  a  double  cantilever  beam  in  the  presence  of  bridging  mechanisms. 
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Figure  8:  Dimensionless  critical  load  versus  normalized  crack  length  in  orthotropic  DCB  specimens:  a) 
No  bridging  mechanisms,  p30has  /^icEl  =  0.0;  b)  pi0h0>  / =  0.1;  c)  pxh05  /Jg =  0.3. 

(curves  (a):  A  =  1  and p=  1;  curves  (b)  A  =  0. 1  and p—  3;  curves  (c)  A  =  0.05  and p=  5) 
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ABSTRACT 

The  dynamic  response  of  laminated  beam-type  structures  presenting  delaminations  is  characterized  by 
stiffness  degradation,  natural  frequency  shifts  and,  in  special  cases,  the  occurrence  of  opening  modes 
during  vibration.  The  paper  investigates  how  bridging  mechanisms  acting  between  the  surfaces  of  the 
delaminations  may  affect  these  behaviors.  Such  mechanisms  could  be  developed  by  a  through¬ 
thickness  reinforcement  (stitching,  z-pins)  applied  to  the  laminate  in  order  to  improve  its  resistance  and 
tolerance  against  delamination  fracture.  A  model  based  on  the  theory  of  bending  of  laminated  plates  is 
formulated  where  the  bridging  mechanisms  are  modeled  as  a  uniform  distribution  of  linear  elastic 
springs  that  oppose  longitudinal  and  transversal  relative  displacements  between  the  surfaces  of  the 
delamination.  The  model  is  applied  to  predict  natural  frequencies  and  modal  shapes  of  a  delaminated 
cantilever  beam  with  a  rectangular  cross  section.  A  variety  of  solutions  and  transitions  in  the  dynamic 
response  is  found  on  varying  the  stiffness  of  the  ligaments.  A  transition  from  damaged  (unbridged)  to 
intact  beam  values  is  predicted  upon  increasing  the  stiffness  of  both  longitudinal  and  transversal 
springs.  An  application  of  the  model  to  a  stitched  carbon-epoxy  laminate  for  the  aeronautical  industry 
shows  that  low  percentages  of  through-thickness  reinforcement  can  substantially  improve  the  dynamic 
response  of  delaminated  structures. 

Keywords:  Delamination,  dynamics,  natural  frequencies  of  vibration,  beam  theory,  composite 
laminates,  through-thickness  reinforcement. 


1.  INTRODUCTION 

Fiber  reinforced  polymer  matrix  laminates  are  characterized  by  high  strength  and  stiffness  to  weight 
ratios  that  make  them  particularly  suited  to  develop  load-bearing  components  for  a  variety  of 
applications  of  the  aerospace,  naval  and  aeronautical  industries.  The  utilization  of  these  materials. 
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however,  has  been  limited  by  their  high  sensitivity  to  interlaminar  flaws  that,  in  the  presence  of 
interlaminar  stresses,  can  easily  propagate  and  lead  to  catastrophic  failures  or  unacceptable  loss  of 
stiffness.  Interlaminar  flaws  may  originate  during  fabrication,  mainly  due  to  air  trapped  in  the  resin  and 
also  to  imperfect  bonding,  matrix  cracks  or  broken  fibers;  they  may  also  be  service-induced,  as  because 
of  impact  with  foreign  objects  and  fatigue  loading. 

Besides  being  a  major  cause  of  stiffness  and  strength  degradation,  delamination  also  causes 
changes  in  the  dynamic  response  of  composite  laminates.  Extensive  experimental  and  theoretical 
investigations  have  shown  that  delaminations  cause  shifts  of  the  natural  frequencies  of  vibration  along 
with  changes  in  the  modal  shapes.  Openings  of  the  delaminated  surfaces  during  vibration,  the  so-called 
“opening  modes”,  have  also  been  experimentally  detected  and  theoretically  reproduced  [1-3]. 

In  this  paper,  a  model  based  on  a  particularization  of  the  first  order  shear  deformation  theory  for 
anisotropic  laminated  plates  in  cylindrical  bending  [4-6]  has  been  formulated  in  order  to  study  the 
dynamic  response  of  delaminated  composite  plates  in  the  presence  of  bridging  mechanisms  acting 
between  the  crack  surfaces.  Following  the  approach  described  in  [3],  a  through-width  delaminated  plate 
has  been  modeled  by  using  four  elementary  plates  connected  at  the  delamination  edges.  The  bridging 
mechanisms  are  represented  by  a  uniform  distribution  of  springs  that  oppose  longitudinal  and 
transversal  relative  crack  displacements  [7].  In  this  paper  only  linear  elastic  springs,  exerting 
symmetric  restoring  forces  during  vibration,  have  been  considered.  This  assumption  leads  to  analytical 
or  semi-analytical  solutions  of  the  problem  allowing  a  complete  description  of  the  variety  of  modal 
shapes  and  associated  natural  frequencies  that  arise  due  to  the  presence  of  bridging  ligaments  of 
different  stiffness.  Nonlinear  bridging  mechanisms  and  the  problem  of  nonlinear  contact  during 
vibration  will  be  studied  elsewhere. 

Bridging  mechanisms  acting  between  the  crack  surfaces  could  represent  the  action  developed  by  a 
through-thickness  reinforcement,  such  as  stitching  or  z-pins.  The  efficacy  of  the  technology  of  through¬ 
thickness  reinforcement  in  enhancing  delamination  resistance  and  damage  and  impact  tolerance  of 
composite  laminates  has  been  abundantly  documented  in  the  last  twenty  years  [8-11].  The  through¬ 
thickness  reinforcements  create  large  zone  of  bridging  along  the  wake  of  the  delamination  that  shield 
the  crack  tip  from  the  applied  loads  thus  reducing  the  driving  force  for  crack  propagation  or  even 
suppressing  crack  growth  [7,12-15].  A  through-thickness  reinforcement  can  reduce  the  loss  of  stiffness 
due  to  delamination,  modify  modes  of  failure  and  prevent  catastrophic  collapses  of  the  structural 
components. 

It  is  expected  that  a  through-thickness  reinforcement  will  also  reduce  the  delamination-induced 
degradation  of  the  dynamic  properties  of  composite  structures  by  opposing  the  relative  displacements 
between  the  delamination  surfaces  during  vibration.  This  behavior  will  be  investigated  using  the 
proposed  theoretical  model. 

The  analysis  of  the  degradation  of  the  dynamic  properties  of  delaminated  beam-type  structures  in 
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the  presence  of  a  through-thickness  reinforcement  has  special  relevance.  It  is  well  known  that  the 
influence  of  short  delaminations  on  the  first  natural  frequencies  of  vibration  is  limited  and  that  only 
considerably  long  delaminations  degrade  the  first  frequencies  to  a  certain  extent  [1-3].  Such  long 
delaminations,  however,  would  not  be  acceptable  in  traditional  laminate  structures.  Through-thickness 
reinforced  laminates,  on  the  other  hand,  can  be  used  within  a  damage  tolerant  design  of  composite 
structures.  These  damage  tolerant  structures  would  be  able  to  tolerate  delaminations  of  considerable 
length  due  to  the  action  developed  by  the  through  thickness  reinforcement.  It  is  then  important  to 
understand  and  quantify  the  effects  these  long  delaminations  might  have  on  the  dynamic  response  and 
investigate  the  influence  of  the  through-thickness  reinforcement. 


2.  MODELING  AND  TESTING  OF  DELAMINATED  BEAMS 

Free  vibrations  of  delaminated  composite  beams  were  first  studied  by  Ramkumar  et  al.  in  1979  [16]. 
Their  analytical  model,  based  on  four  Timoshenko’s  beams  connected  at  the  delamination  edges  (Fig. 
1),  neglected  the  coupling  between  flexural  and  axial  vibrations  and  therefore  led  to  an  overestimation 
of  the  drop  in  the  natural  frequencies  caused  by  delamination  and  observed  in  experimental  tests.  Wang 
et  al.  [1]  improved  the  analytical  model  by  including  the  coupling  effect  and  obtained  a  better 
agreement  with  experimental  frequency  measurements.  In  these  first  models,  also  known  as  free 
models,  the  two  sub-laminates  above  and  below  the  delamination  plane  were  allowed  to  vibrate 
independently.  For  special  geometrical  configurations  (very  long  cracks  close  to  the  surface  of  the 
beam)  this  assumption  could  lead  to  overlapping  of  the  sublaminates  during  portions  of  the  vibration 
cycle. 

To  avoid  the  overmentioned  incompatibility,  Mujumdar  and  Suryanarayan  [17]  imposed  the  two 
delaminated  parts  to  have  the  same  transversal  displacement  during  vibration.  This  so-called 
constrained  model  was  later  adopted  by  Tracy  and  Pardoen  [18]  to  study  the  free  vibrations  of  simply 
supported  isotropic  beams  and  showed  good  agreement  with  finite  element  results  and  experimental 
frequency  measurements. 

The  problem  of  the  interaction  between  sublaminates  during  free  and  forced  vibrations  of 
composite  beams  in  the  presence  of  one  or  more  delamination  cracks  has  been  investigated  in  [3,19- 
23].  Luo  and  Hanagud  [3]  studied  the  problem  of  a  beam  with  a  single  delamination  applying 
Timoshenko  beam  theory.  They  modelled  the  contact  forces  that  arise  during  the  compression  cycles 
using  a  uniform  distribution  of  transversal  springs.  A  piece-wise  linear  spring  stiffness  was  defined  to 
allow  transversal  relative  displacements  up  to  a  pre-set  amount  and  simulate  the  transversal  elasticity  of 
the  two  sublaminates.  Further  overlap  was  impeded  by  an  infinite  spring  stiffness.  The  springs  opposed 
no  resistance  to  opening  displacements.  Different  frequencies  and  modal  shapes  were  then  defined 
corresponding  to  the  different  spring  stiffness  and  the  results  combined  in  the  solution  of  a  forced 
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vibration  problem.  Lee  in  [22]  applied  a  higher-order  plate  theory  within  a  layerwise  approach 
[24,25,26]  to  numerically  investigate  the  free  vibrations  of  delaminated  plates. 

Other  recent  studies  have  dealt  with  the  problem  of  free  vibrations  with  respect  to  post-buckled 
states  [27-29],  and  the  problem  of  the  influence  of  delamination  cracks  on  the  flutter  characteristics  of 
laminated  plates  [30].  Flutter  is  a  self-excited  oscillation  caused  by  the  aerodynamic  forces  generated 
by  the  motion  of  a  structure  immersed  in  a  wind  flow  -  once  a  critical  velocity  has  been  reached,  the 
structure  starts  to  vibrate  violently.  In  [30]  the  presence  of  delaminations  has  been  shown  to  alter  the 
conditions  characterizing  the  onset  of  flutter  in  composite  plates  and  reduce  the  flutter  boundary. 


3.  ASSUMPTIONS  AND  BASIC  FORMULATION  OF  THE  MODEL 


Figure  1  shows  a  laminated  plate  of  length  L,  thickness  t  and  width  w  (w»/)  with  a  through-width 
delamination  of  length  Ldei  parallel  to  the  x-y  plane,  with  x  the  longitudinal  geometrical  axis  of  the 
plate.  The  delamination  identifies  four  plates,  for  i=  1,  ...4,  of  length  L„  thickness  t„  cross  sectional 
area  Si  and  central  moment  of  inertia  /j.  The  distance  of  the  plate  centroidal  axes  from  the  delamination 
plane  is  indicated  with  h\.  The  position  of  the  mid-point  of  the  delamination  in  the  plane  x-z  is 
identified  by  Lc  =  Li+  L& i/2  and  h. 

The  laminate  consists  of  n  linearly  elastic  and  transversally  isotropic  laminae  with  material  density 
p .  The  principal  material  axes  of  the  laminae  are  arbitrarily  oriented  and  the  layup  is  such  to  ensure 
plain  strain  conditions  parallel  to  the  plane  x-z.  Only  small  deformations  and  small  displacements  are 
considered.  Based  on  these  assumptions,  the  plate  vibrates  in  cylindrical  bending  and  the  governing 
equations  coincide  with  those  of  a  beam  characterized  by  a  reduced  flexural  stiffness.  In  the  following, 
the  terms  beam  and  plate  will  be  used  interchangeably. 

Linear  bridging  mechanisms  acting  along  the  surfaces  of  the  delamination  are  represented  in  the 
model  by  means  of  uniformly  distributed  linear  elastic  springs  that  exert  symmetric  restoring  forces  to 
both  transversal  and  longitudinal  relative  displacements  between  the  delamination  surfaces.  The 
stiffness  in  the  transverse  and  longitudinal  directions  are  indicated  with  kv  and  kt,  respectively. 

The  generalized  displacements  of  the  axis  of  /th  beam  are  the  transverse  displacement,  w>i,  the  axial 
displacement,  u„  and  the  bending  rotation,  %  the  shear  deformation  y{  is  related  to  wj  and  %  through 
the  compatibility  relation: 


Y\  =  <Pi  + 


dwt 

dx 


(1) 


The  stress  resultants  per  unit  width  are  the  normal  force,  N„  shear  force,  Q„  and  bending  moment,  M-,t 
for  i  =  1,..,  4  (Fig.l). 

The  constitutive  equations  for  the  ith  beam  are: 
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Ni=A^+BiKxi  ,  (2) 

Mi  =  +  A^xi  .  (3) 

Q-^G-Jy  ,  (4) 

where  =  uifX  and  Kxi  =  ^JjJt  are  the  axial  deformation  and  bending  curvature  and  the  comma 
indicates  partial  differentiation;  the  laminate  stiffness  coefficients  A„  B„  D\  and  Gi  are  defined  using 
lamination  theory  and  the  assumptions  in  [6]. 

The  equations  of  motion  for  the  intact  parts  of  die  beam  are: 

^i»x  =  /^iWi>tt  > 

Qi,x  =  yaS’iwi,tt  ,  i=l,4  (5a,b,c) 

Mi,x-Qi=pIi<Pi,n  , 

where  the  terms  on  the  right  hand  side  account  for  the  inertia  of  the  plate  and  the  symbol  t  has  been 
used  for  the  time  variable.  For  the  upper  and  lower  sublaminates,  the  equations  of  motion  are: 

^2»x+^/(W3  —  u2  ~<P2^l2  ~  ^3^3  )=  P^2M2’tt  » 

22»x+^v(w3  -w2)  =  P52w2>tt  »  (6a,b,c) 

M29X~Q2+kth2(u3—U2—^>2^t2~P3^h)~P^2<P2>n  * 

^3»x-^(«3-"2-^2^2-?,3A3)=/jS'3«3»tt  » 

03»x-^v(w3-w2)=/^3w3»tt  »  (7a,b,C> 

^3>x-03  +kth i(u3  -U2  -^2^2  ~(P^rll)~  • 

Equations  (5-7)  can  be  expressed  in  terms  of  displacement  components  by  exploiting  equations  (2-4) 
and  the  compatibility  relations. 

A  harmonic  type  solution  is  assumed  for  the  free  vibration  problem.  The  solution  is  characterized 
by  a  single  angular  frequency  (O  and  the  displacement  field  assumes  the  following  form: 
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wi(x,f)  =  wi(x)ey‘0'  Wj(x,r)=  Wi(x)e-"“  <p[{x,t)=(pXx)eio*  i=l, -4.  (8a,b,c) 

Substitution  of  Eqs.  (8a,b,c)  into  the  equations  of  motion  leads  to  a  system  of  twelve  ordinary 
differential  equations.  To  determine  the  twelve  displacement  amplitudes,  tq(x),  ^(x),  <pt (x),  and  the 
angular  frequency,  a ),  the  eigenvalue  boundary  value  problem  is  completed  by  continuity  and  boundary 
conditions.  Such  equations  will  be  explicitly  stated  in  a  dimensionless  form  in  the  following  section. 

4.  DIMENSIONLESS  FORMULATION 

The  thickness  of  the  intact  beam,  t,  its  extensional  stiffness  coefficient,  A\,  and  the  material  density,  p, 
are  chosen  as  fundamental  quantities  to  define  the  following  dimensionless  groups  for  the  geometric 
variables,  the  generalized  displacements  and  the  elastic  stiffness  coefficients: 

t2=l-t3  ,  it-,  =~~> 

(9a,b,c) 


with  i  =  1,  ...4.  The  dimensionless  stiffness  parameters  for  the  transversal  and  longitudinal  springs  are: 

k  t2  k  t2 

a  =  5d_  and  p  =  ,  (9d,e) 

A\  Ax 

and  the  dimensionless  eigenvalue  is: 

(9f) 

At 

The  twelve  governing  equations  in  the  thirteen  unknowns  u, ,  w, ,  <pt  (i  =  1....4)  and  y/  are  expressed 
in  a  dimensionless  form  as  follows: 

A 

i=1>4  *  (10a,b) 
Gi 
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-  I  1  1  a  ¥\V  i  I  -  Bi  «  a  .  i  . 

w.x,ms+¥  7r+—K  wltfe+7-  — y-1  wt-j-ux^ 0  1-1,4  , 

V  Gj  \2DiJ  DA12G;  J  Di 


(Ua,b) 


V  «  .  i  . 

pi)5  =  -4-w i-w,,#  1=1,4  , 

C T, 


(12a,b) 


*  ^2  a  B 2  «  -®2^2  *  -®2  / a 

"2^+^^«2^-l-W2^-^l-7TW2^-«TX-K> 


P  A  A  ^2^2  A  ,*  A  ^3  /j  A  r  A  r^2  ^3  /a  A  \  A 

-y  ui,4-u2,4+\^-^-w2  +  h2w2,g+y/-^f-w3  +h3w3,44+a\  -jr—4-.  K  -w2)  =  0 


G2  g3 


,  (13) 


G2  12Z>2 


•«WHrf2[4-  +  +^(*3.«-W2.«)+-y|yf-l  p3  -*2)- 


a 

^  A  O 

'2V 

4 

[l2G2 

A  A  A  A  A  A  \ 

A  A  A  .A  A  A  •‘a  ^9  /ll  1/  A  A  \  A 

W3^-W2,|+^-^W2+^2w2^+^-^-W3+;*3w3^+a  J—J-h  “W2)  =0 


G2  G3 


.  (14) 


^  A  a  ,  A  A  \  A 

<P2>E  =  — rW2- S-lW3-W2)-W2,»  , 


A  ^3  A  A  -®3^3  A  -63  /  A  A  \ 

«3  .wr+i^  -y  "3  *  -y  w3»«k  ~^~Ta~W3  ’«  +aT7-^3  ’#  _W2,« 

ilj  iT-3  AJU3  il^Uj 

A  A  A  ^  /a  A  \ 

P  ~  A  h^tj  A  ,  f  A  ,  ^3^3  A  ,  /“  A  .  h2  h3  /a  A  \  A 

-y  «3»^-“2^+^-^W2+A2w2.||+^-y-w3+A3W3»^+«  \W2~W2)  =° 

•"3  @2  G3  \G2  ^3  y 


,  (16) 


G3  12£>3 


4 

f  a3 

'3V 

a 

v12G3 

a  /a  a  ^ 

(  A3 

/3V 

7Tvw3»^“w2^J 

G3 

1-— 

4 

,12G3 

/^3  A  A  ^2^2  *  f  A  ^3^3  A  r  A  ^2  ^3  |/a  A  \  « 

-y  K3^-«2.f+^-y-w2+/»2M'2»«+^yi-M'3+/»3w3»ff+«  -y-~y  (w3 -^2)  =° 


g2  g3 


g3  g3 

The  kinematic  and  static  continuity  conditions  at  the  delamination  edge  are: 


7 


Submitted  to  Composite  Structures,  October  2002.  Accepted,  December  2002. 


A  A 

w,  =w2. 

A  A 

W]  =w3> 

(19a, b) 

9S 

& 

II 

£ 

£ 

II 

(20a, b) 

«2  =w,  +  (a,  -h2)p ,, 

«3  =«,+(a,+A3)p1, 

(21a, b) 

M,=  M2  +  M3-  (h2  -h^N.+fa+h^N,,  0,  =  02  +  03,  Nl=N2+N2  .  (22a, b,c) 

Analogous  relations  are  written  for  the  other  delamination  edge  at  £=  (Lj  +  Ldel  )/t .  Equations  (21a,b 

and  22a)  refer  to  all  cases  where  the  centroidal  axis  of  the  beam  1  is  above  the  plane  of  the 
delamination.  Cases  where  the  centroidal  axis  of  the  beam  1  is  below  the  delamination  plane  are 

A  A 

described  by  the  same  equation  substituting  A,  (or  A,)  with  -  A,  (or  -  A,). 

For  a  cantilever  beam  clampled  at  £  =  0  the  boundary  conditions  in  the  absence  of  applied  forces 
are: 

«,  =  0,  w,  =0,  <px  =0  ,  at£  =  0  »  (23a,b,c) 

N4=0,  Ma  =0  ,  04  =0  ,  at  Z  =  L  .  (24a, b,c) 

The  system  of  equations  (10-18)  and  the  continuity  and  boundary  conditions  (19-24)  constitute  a 
two-point  boundary  eigenvalue  problem  with  internal  conditions  in  terms  of  the  generalized 
displacements  wt,  <pt  and  the  eigenvalue  y .  The  problem  generally  requires  a  numerical  integration 
such  as  that  presented  in  Appendix  I. 

The  proposed  model  has  been  validated  by  referring  to  the  asymptotic  frequency  values  obtained  by 
letting  the  spring  stiffness  vanish  (thus  simulating  the  absence  of  reinforcement)  and  indefinitely  grow 
(absence  of  delamination).  For  a  vanishing  spring  stiffness,  the  experimental  tests  of  Shen  and  Grady 
[2],  on  delaminated  cantilever  beams,  have  been  reproduced  [31,32].  Luo  and  Hanagud  [3]  used  the 
same  tests  to  validate  their  beam  model.  For  very  high  spring  stiffness,  the  frequency  values  predicted 
by  the  model  approach  the  values  corresponding  to  an  intact  beam  [33]. 


5.  ANALYTICAL  SOLUTION 

A  closed-form  solution  of  the  boundary  eigenvalue  problem  presented  in  the  previous  section  is 
obtained  under  the  following  additional  assumptions: 
a)  mid-plane  delamination:  t2=t3=t/2; 
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b)  constrained  model:  w2  =  w3  =  w ; 

c)  homogeneous  material; 

d)  negligible  shear  deformations:  y,=  0  for  i  =  1,...4. 

As  it  will  be  shown  later,  the  assumption  (b)  is  satisfied  in  all  cases  but  for  very  long  delaminations 
close  to  the  surface  of  the  beam  when  opening  modes  may  develop  [1-3].  In  all  other  cases  the 
transversal  springs  are  not  activated  during  vibration  and  the  solution  therefore  does  not  depend  on  the 
normalized  stiffness  coefficient  a  (Eq.  9d).  The  assumption  (c)  is  satisfied  by  all  isotropic  and 
orthotropic  homogeneous  materials  and  by  unidirectional  lay-ups  made  of  the  same  orthotropic 
laminae.  Furthermore,  laminates  with  many  layers  and  no  bending/extensional  coupling,  such  as  quasi¬ 
isotropic  laminates  typically  used  for  aeronautical  applications  and  cross-plies  [0/90]n,  belong  to  this 
category.  The  assumption  leads  to  Bt  =0  (i  =  1....4)  and  A,  =Elt,  with  the  reduced  longitudinal 
Young’s  modulus. 

The  problem  reduces  to  the  determination,  besides  the  eigenvalue  y/ ,  of  the  seven  functions  w,(|), 
w4(|)  and  w(£),  for  i  =  representing  the  dimensionless  displacement  amplitudes. 

Equations  (10,11,13,16)  simplify  as  follows: 


u^g+ydii  =0  ,  i  =  1,4  , 

(25a,b) 

W;  ,m  +yrii> j  ,R  -\2y/w{  =  0  ,  i  =  1 ,4  , 

(26a, b) 

(  w,£  'l 

w3-u2+  — 1  =  0  , 

(27) 

«3-«2+^iJ  =  0  •  (28) 

Moreover,  by  introducing  die  relative  longitudinal  displacement  ut  =  u3  -  u2 ,  the  last  two  equations 
can  be  combined  to  give: 

ut,g+(yr-4fi)ut  +  2fiw,s  =  0  .  (29) 

Finally,  Eqs.  (14)  and  (17)  reduce  to  the  following  expression: 

-(l  2fi  -  y/)w,  s  -48^  =  24 p  m,  .  (30) 

Equations  (25, 26, 29  and  30)  constitute  a  system  of  six  ordinary  differential  equations  in  the  unknowns 
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y/,  «i(|),  u4(£),  w,(|),  w4(4),  ut(t)  and  w(£).  The  kinematic  continuity  conditions  (19-20)  at  the 
first  delamination  edge  4i  become: 


Wi=W,  = 


-  ™i>{ 


«3  =  «i 


wi»f 


(31) 


Similar  conditions  are  posed  at  the  other  delamination  edge.  The  static  continuity  conditions  (22a, b,c) 
and  boundary  conditions  (23)  and  (24)  remain  unchanged. 

The  solution  procedure  of  the  system  of  equations  is  outlined  in  the  Appendix  II.  The  general 
solution  is: 

«i  (£)=<:,  sin  (<?,£)  , 
w4(£)=C2cos[<5,  (£-&)]  , 

w, (£)  =  C3  sin  (tj,4)+  C4  cos  (r)x4)+  C5  sinh  (<?,£)+  C6  cosh  (0X4)  , 

w4(4)=C7 sin (tjx4)+  Cg cos (r}x4)+  C9 sinh {Ox4)+  C10 cosh (<?,£)  , 

«,(£)  =  Cnfx  sin  {ax4)-Cnfx  cos  (a,£)+  C13g2  cosh  (a2<*)- CI4g2  sinh 
+  C15g3  cosh  (a34)  +  C16g3  sinh  (a3£) 

w{4)  =  Cn  cos  {ax4)+  Cx  2  sin  (a,£)+  C13  cosh  (a2£)+  C14  sinh  (a2<f)+ 

+  C15  cosh  («3^)+  C16  sinh  (a3£) 

— — u2(4)  =  CnPx  sin  (a,#)-  C12?,cos  («,£)+  CX3p2cos  (a2£)-  CX4q2 sin  (a2^)+ 

+  C,5/73cosh(a3#)+C16g3sinh(a3^)+C17cos(^)+C1gSin(i/^) 

where  the  parameters  , rjj ,3j,ai,  fx, px , qx  and C,  are  defined  in  the  Appendix  II.  They  depend  on 
geometrical  parameters,  material  properties  and  the  natural  angular  frequency  0).  The  continuity  and 
boundary  conditions  (22-24,31)  define  a  system  of  eighteen  algebraic  equations  in  the  eighteen 
integration  constants  C„  i  =  1,  ..18. 

The  natural  frequencies  and  corresponding  modal  shapes  are  then  calculated  by  imposing  the 
vanishing  of  the  determinant  of  the  system  coefficients.  It  is  worth  noting  that  for  very  low  and  very 
large  values  of  the  normalized  longitudinal  stiffness  of  the  ligaments,  /? ,  a  high  precision  is  required  in 
the  numerical  solution. 


(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


10 


Submitted  to  Composite  Structures,  October  2002.  Accepted,  December  2002. 


6.  NUMERICAL  SOLUTION:  CONSTRAINED  MODEL 

If  any  of  the  assumptions  stated  at  the  beginning  of  the  previous  section  are  removed,  a  closed  form 
solution  of  the  two-point  boundary  eigenvalue  problem  is  no  longer  attainable.  It  is  still  possible 
however  to  proceed  with  a  numerical  integration  using  the  method  outlined  in  the  Appendix  I,  which  is 
based  on  a  relaxation  technique  [34]. 

The  model  described  in  Sections  3  and  4  has  been  used  to  investigate  the  role  played  by  the 
bridging  mechanisms  in  reducing  the  delamination-induced  degradation  of  the  dynamic  properties  of 
delaminated  cantilever  beams.  In  this  section,  results  obtained  under  the  assumptions  of  constrained 
model,  w2  =  w3  =  w,  homogeneous  material  and  negligible  shear  deformations,  will  be  presented  in 
order  to  investigate  the  influence  of  the  longitudinal  springs.  In  the  next  section,  a  special  geometry 
leading  to  opening  modes  will  be  examined  to  investigate  the  influence  of  the  transversal  springs. 

The  results  are  presented  in  terms  of  normalized  natural  frequency  of  vibration  ,  with^o  the 
natural  frequency  of  the  corresponding  intact  beam,  dimensionless  longitudinal  stiffness  of  the 
ligaments  p  =  ktt/Ex,  normalized  length  of  the  delamination,  Liel/L,  normalized  length  of  the  beam, 
L/t,  and  normalized  position  of  the  mid-point  of  the  delamination,  Lc/t .  The  parameter  ft  has  been 
varied  between  1 0‘7  to  1.0. 

It  is  well  known  that  the  first  natural  frequency  of  vibration  is  degraded  only  by  very  long 
delaminations.  The  diagram  of  Fig.  2,  which  refers  to  LM/L=  0.50,  L/t  — 50  and  Lc/L  =  0.5  and 
shows  the  first  natural  frequency  versus  the  dimensionless  stiffness  ft  for  different  thicknesswise 
positions  of  the  delamination,  confirms  this  result.  For  very  small  values  of  P  all  curves  tend  to  the 
normalized  frequencies  of  the  delaminated  beam  in  the  absence  of  through-thickness  reinforcement  [3]. 
As  expected,  the  curves  show  a  transition  from  damaged  and  unbridged  values  to  the  values 
corresponding  to  an  intact  beam  (f=fo)  upon  increasing  the  elastic  spring  stiffness  p .  The  maximum 
degradation  caused  by  the  delamination  corresponds  to  the  case  t2/t  =  0.5  -  for  higher  values  of  t2 /t 
the  dynamic  properties  of  the  delaminated  beam  are  progressively  more  affected  by  the  behavior  of  the 
thicker  sublaminate. 

The  chart  in  Fig.  3  illustrates  how  the  ///„  ratio  corresponding  to  the  first  natural  frequency  varies 
as  a  function  of  P  for  different  delamination  locations  along  the  beam  length.  The  chart  refers  to 
Am  /L  =  0.50 ,  t2/t  =  0.5  and  L/t  =  50  (low  values  of  LjL  indicate  delaminations  close  to  the 
clamped  end).  The  frequency  reduction  is  stronger  for  delaminations  close  to  the  clamped  end. 

Figure  4  refers  to  the  first  three  natural  frequency  ratios  plotted  as  a  function  of  P  for  a  beam  with 
L/t  =  50,  LjL  =  0.5 ,  Ldel  /L  =  0.25  and  t2  /t  =  0.5 .  It  can  be  noted  that  the  reduction  of  the  third 
natural  frequency  is  considerably  higher  than  those  corresponding  to  the  first  two  modes.  Likewise,  the 
effectiveness  of  the  through-thickness  reinforcements  is  correspondingly  more  evident,  taking  the  value 
of  the  third  natural  frequency  from  0.843/o  tofo  on  increasing  p. 
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It  is  important  to  note  that  for  all  results  presented  in  this  section,  apart  from  the  curve  of  Fig.  2 
corresponding  to  t2/t  =  0.9  (see  next  section),  the  approximation  of  constrained  model  coincides  with 
the  rigorous  solution  and  the  sublaminates  above  and  below  the  delamination  plane  vibrate  with  the 
same  transversal  displacement. 

7.  NUMERICAL  SOLUTION:  GENERAL  CASE 

In  the  results  presented  in  this  section  the  assumption  of  constrained  model  has  been  relaxed  allowing 
the  two  sublaminates  above  and  below  the  delamination  plane  to  vibrate  independently.  The 
dimensionless  transversal  spring  stiffness  a  has  been  varied  between  two  limit  values,  as  and  ac .  The 
two  values  have  been  chosen  so  that  when  a  tends  to  as  the  frequency  ratio  approaches  the  value 
corresponding  to  a  =  0,  a  condition  referred  to  as  shear-bridged  model  that  coincides  with  the  free 
model  in  the  special  case  /?  =  0;  when  a  tends  to  ac  the  frequency  ratio  approaches  the  value 
obtained  through  the  constrained  model  and  further  increases  in  the  value  of  a  do  not  modify  the 
solution. 

Following  the  experimental  and  theoretical  observations  in  [2,3],  a  slender  beam,  L/t  =  150,  with  a 
long  delamination,  LitXjL  -  0.5 ,  close  to  the  surface,  t2jt  =  0.9 ,  has  been  examined  in  order  to  be  able 
to  observe  the  phenomenon  of  crack  opening  during  vibration.  The  delamination  is  centrally  located, 
Lc/L  =  0.5.  As  in  the  previous  section,  results  are  presented  for  a  homogeneous  material  and 
negligible  shear  deformations. 

In  the  presence  of  opening  modes,  the  closure  forces  exerted  by  the  transversal  springs  when  the 
delamination  opens  represent  the  action  of  the  bridging  mechanisms,  whereas  the  opening  forces  of  the 
same  springs  simulate  the  elastic  contact  forces  due  to  the  closure  of  the  delamination.  Therefore, 
allowing  relative  transversal  displacements  of  the  sublaminates  causing  springs  compression  does  not 
necessarily  imply  overlapping  of  the  sublaminates  axes,  as  already  observed  in  Mujumdar  et  al.  [16],  as 
long  as  the  vibration  amplitudes  remain  sufficiently  small.  In  this  paper  the  spring  stiffness  has  been 
assumed  constant  in  the  tensile  and  compressive  regimes.  This  assumption  simplifies  the  solution  of 
the  problem  and  is  not  expected  to  substantially  affect  the  conclusions. 

The  influence  of  the  spring  stiffness  parameters  on  frequency  ratios  and  modal  shapes  of  a 
cantilever  beam  are  described  in  Sections  7.1  and  7.2,  respectively.  Due  to  the  slenderness  of  the  beam, 
the  solutions  presented  in  this  section  will  be  quantitatively  correct  only  in  a  regime  of  vibration 
characterized  by  small  amplitudes. 

7.1  Influence  of  spring  stiffness  parameters  on  frequency  ratios 

The  frequency  ratio  ///0  corresponding  to  the  first  natural  frequency  of  vibration  is  plotted  in  Fig.  5 
versus  the  longitudinal  stiffness  f  for  as  =KT10,  ac  =10-8  (corresponding  to  the  shear-bridged  and 
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constrained  models,  respectively)  and  an  intermediate  value  a  =  10”9.  It  can  be  noted  that  the  solution 
is  only  slightly  affected  by  a  and  the  intact  beam  frequency  ratio  is  always  approached  by  adequately 
increasing  /? .  The  negligible  influence  of  a  can  be  explained  by  noting  that,  even  in  the  absence  of 
bridging  mechanisms,  the  first  modal  shape  of  the  delaminated  beam  is  characterized  by  a  small 
opening  of  the  delamination,  thus  slightly  differing  from  the  constrained  model  condition  (see  Fig.  1 1). 

For  the  second  natural  vibration  mode,  the  frequency  reduction  induced  by  the  delamination  on  the 
cantilever  beam  in  the  absence  of  bridging  mechanisms  is  59.2%  and  the  limits  on  a  are  as  =  10-11  and 
ac  =  10-7.  Figure  6  shows  a  plot  of  the  frequency  ratio  versus  p  for  a  =  as,  a  =  ac  and  three 
intermediate  values.  In  this  case  the  intact  beam  frequency  ratio  is  approached  only  for  a>ac.  As  it 
will  be  shown  in  Section  7.2  (Fig.  12),  the  modal  shape  of  the  delaminated  beam  in  the  absence  of 
bridging  mechanisms  strongly  differs  from  the  constrained  model  approximation,  being  characterized 
by  a  delamination  opening  whose  amplitude  is  considerably  higher  than  that  defining  the  vibrational 
motion  of  the  rest  of  the  beam;  therefore,  for  a  <ac  the  condition  // /„  =1  cannot  be  approached. 

The  trend  described  for  the  second  natural  frequency  occurs  for  the  third  natural  mode  as  well.  In 
this  case,  the  frequency  reduction  induced  by  the  delamination  on  the  beam  in  the  absence  of  bridging 
mechanisms  amounts  to  63.4%  and  the  limits  on  a  are  as  =10-9  and  ac  =10'*.  The  plot  of  the 
frequency  ratio  versus  (5  for  a  =  as,  a  =  ac  and  three  intermediate  values  is  shown  in  Fig.  7.  In  this 
case  the  frequency  ratio  of  an  intact  beam  is  no  longer  approached  by  simply  raising  a  even  for 
a>ac.  The  chart  also  highlights  a  limited  effect  of  P  in  all  cases  where  a  <ac. 

Figures  8-10  show  the  dependence  of  the  frequency  ratio  on  the  transversal  spring  stiffness 
parameter  a .  The  curves  in  Fig.  8  confirm  the  negligible  influence  of  a  on  the  first  frequency  ratio. 
The  chart  of  Fig.  9,  referring  to  the  second  natural  frequency,  clearly  shows  that  the  frequency 
reduction  induced  by  the  delamination  is  completely  recovered  only  if  a  >  ac  =  10~7.  For  decreasing 
values  of  a ,  the  curves  show  a  growing  spread  in  the  vertical  direction,  which  confirms  the  higher 
effectiveness  of  the  longitudinal  springs.  Similarly,  the  curves  in  Fig.  10  show  that,  for  the  third  natural 
frequency,  the  degradation  induced  by  the  delamination  is  almost  recovered  if  a  >ac  =10"®  and 
P  >  10"4 .  In  this  case,  however,  the  curves  are  practically  unaffected  by  y 6  for  low  values  of  a . 

7.2  Influence  of  spring  stiffness  parameters  on  modal  shapes 

The  two  charts  in  Figure  11  show  the  first  modal  shape  of  the  cantilever  delaminated  beam 
corresponding  to  the  two  limiting  solutions  of  constrained  model  with  /?=  0  and  free  model  (a  =  0  and 
P=  0).  In  these  graphs  and  those  that  follow,  the  modal  shapes  have  been  normalized  assuming  a  unit 
deflection  at  the  free  end.  The  horizontal  and  vertical  axes  define  the  normalized  position  along  the 
beam,  E,=x/t,  and  the  normalized  deflection,  w  =  w(x)/w(L).  For  a  =  1CT10  the  presence  of  the 
delamination  causes  an  almost  undetectable  opening  mode  that  disappears,  as  expected,  for 
a  =  ac  - 10-8 .  However,  since  P  =  0,  the  normalized  frequency  of  vibration  remains  lower  than  one. 
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Increasing  the  value  of  /?  from  0  to  1  leads  to  the  vanishing  of  the  opening  mode  and  a  frequency  ratio 
approaching  unity. 

The  sequence  of  six  charts  of  Figure  12  describes  the  evolution  of  the  second  modal  shape  on 
varying  a  between  asand  ac  with  P  =  0.  For  the  free  model  condition,  corresponding  to  a  =  10_1°, 
the  thinner  sublaminate  shows  a  modal  shape  that  is  almost  symmetric  with  respect  to  the  center  of  the 
delamination,  whereas  the  rest  of  the  beam  vibrates  according  to  the  first  modal  shape  of  an  intact 
cantilever  beam.  By  raising  a ,  the  amplitude  of  the  opening  mode  progressively  decreases  and  the 
rest  of  the  beam  assumes  a  modal  shape  corresponding  to  the  second  natural  frequency  of  an  intact 
beam.  The  effect  of  letting  p  vary  from  zero  to  unit,  for  each  value  of  a ,  is  that  of  reducing  the 
amplitude  of  the  opening  mode  without  introducing  major  alterations  in  the  modal  shape. 

The  evolution  of  the  third  modal  shape  of  the  delaminated  beam  for  increasing  values  of  a  and  for 
P  =  0  is  shown  in  the  nine  charts  of  Figure  13.  Up  to  a  »  2.0  ■  10-8 ,  the  thinner  sublaminate  vibrates 
according  to  the  first  modal  shape  of  an  intact  beam  clamped  at  both  ends  and  is  out  of  phase  with  the 
rest  of  the  beam  that  vibrates  according  to  the  second  modal  shape  of  an  intact  cantilever  beam.  The 
amplitude  of  the  opening  mode  progressively  increases  on  raising  a  and  this  phenomenon  has  to  be 
thought  of  as  the  beginning  part  of  a  more  complex  modification  of  the  modal  shape.  For  a  higher 
then  3.5  •  10~8 ,  the  amplitude  of  the  opening  mode  decreases,  and  two  interesting  modifications  appear 
in  the  modal  shape  of  the  beam:  a)  the  deflected  shape  of  thinner  sublaminate  crosses  the  delamination 
at  an  intermediate  point  along  the  delamination;  b)  the  modal  shape  of  the  rest  of  the  beam  tends 
toward  that  corresponding  to  the  third  natural  mode  of  an  intact  beam.  By  further  increasing  a ,  the 
thinner  sublaminate  progressively  assumes  an  anti-symmetric  modal  shape  across  the  delamination, 
resembling  that  corresponding  to  the  second  natural  mode  of  an  intact  beam  clamped  at  both  ends. 
Finally,  for  a  =  ac  « 10-6,  the  opening  mode  disappears  and  the  beam  vibrates  according  to  a  typical 
third  modal  shape  of  an  intact  beam. 

The  effect  of  P  on  the  third  modal  shape  depends  on  the  value  assumed  for  a .  For  a  <  2.0  - 1 0-8 , 
the  effect  of  raising  p  is  that  of  increasing  the  amplitude  of  the  opening  mode  without  altering  the 
modal  shape  shown  by  the  first  four  charts  of  Figure  13.  On  the  other  hand,  for  a  >  3.5 -10-8,  the 
higher  the  value  of  p ,  the  smaller  the  amplitude  of  the  opening  mode,  even  though  the  constrained 
model  condition  can  be  reached  only  for  a>ac=  10-6.  At  a  «  3.5  •  10“8 ,  the  effect  of  raising  P  is 
illustrated  in  the  three  charts  of  Fig.  14:  the  thinner  sublaminate  vibrates  with  increasing  amplitude 
whereas  the  rest  of  the  beam  is  progressively  constrained  to  modify  its  modal  shape  from  that 
corresponding  to  the  third  natural  mode  of  an  intact  beam  (for  p  =  0)  back  to  the  second  natural  mode 
of  an  intact  beam  (for  /?  =  1). 
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8.  DISCUSSION 

In  order  to  understand  the  variety  of  results  and  transitions  in  the  dynamic  response  observed  in 
Sections  6  and  7,  a  differentiation  must  be  made  between  configurations  characterized  by  the  absence 
of  opening  modes  and  configurations  characterized  by  the  presence  of  opening  modes.  For  the  first  the 
rigorous  solution  of  the  problem  coincides  with  the  constrained  model  approximation. 

Absence  of  opening  modes 

When  the  delamination  is  not  very  long  and  close  to  the  surface  of  the  beam  or  when  the  stiffness 
of  the  transversal  springs  is  higher  than  ac  (i.e„  a  >10”®  -s-lO”6),  the  two  sub-laminates  above  and 
below  the  delamination  plane  vibrate  with  the  same  transversal  displacement.  In  all  these  cases, 
however,  the  presence  of  the  delamination  leads  to  relative  sliding  displacements  between  the  crack 
faces  that,  in  the  absence  of  longitudinal  crack  bridging  mechanisms,  reduce  both  the  flexural  stiffness 
of  the  member  and  the  natural  frequencies  of  vibrations. 

The  different  effect  the  delamination  has  on  the  first,  second  and  third  natural  frequencies  (see  Fig. 
5)  is  explained  by  the  different  sliding  displacements  arising  during  vibration  in  the  corresponding 
modal  shapes.  The  first  and  second  modal  shapes  (see  for  instance  last  charts  of  Fig.  11  and  12)  are 
such  that  very  small  crack  sliding  arises  during  vibration  and  the  beam  behaves  almost  like  an  intact 
beam.  On  the  other  hand,  the  third  modal  shape  is  such  that  large  sliding  displacements  arise  during 
vibration  and  consequently  the  effect  of  the  delamination  on  the  frequency  reduction  is  much  stronger 
(see  for  instance  last  chart  of  Fig.  13). 

The  above  observations  have  been  confirmed  by  numerical  calculations.  The  diagram  ofFig.  15  has 
been  obtained  for  the  special  geometry  of  Section  7  (Liel/L  =  0.5,  LefL  =  0.5,  L/t  =  150  and 
t2  ft  =  0.9)  arid  assuming  fi  =0  in  order  to  describe  how  the  sliding  displacement  between  the  crack 
surfaces  is  controlled  by  the  stiffness  of  the  transversal  springs.  A  parameter  X  has  been  introduced  as  a 
measure  of  the  total  sliding  displacement: 

2= j  -«3  +  <p2h2  +<p3h3\dx.  (39) 

where  the  integrand  defines  the  crack  sliding  at  the  coordinate  x  along  the  delamination.  The  parameter 
X  has  been  calculated  numerically  and  diagrammed  as  a  function  of  a  for  the  first  three  modal  shapes. 
In  the  diagram,  the  absence  of  opening  modes  corresponds  to  values  of  a  higher  than 
orc  =10-8  h-10"6.  The  diagram  shows  that,  for  such  values  of  a ,  X  corresponding  to  the  first  and 
second  modes  of  vibration  is  virtually  zero  while  large  values  of  X  characterize  the  third  mode. 

In  all  configurations  characterized  by  the  absence  of  opening  modes,  the  effect  of  the  longitudinal 
springs  is  always  that  of  restoring  the  conditions  of  an  intact  beam  (see  Figs.  2-4). 
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Presence  of  opening  modes 

In  the  configurations  characterized  by  opening  modes  during  vibration  substantial  reductions  of  all 
natural  frequencies  are  typically  observed.  This  happens  both  when  crack  sliding  displacements  are 
allowed,  i.e.  when  /3  =  0  (see  first  charts  of  Fig.  12-14),  and  when  they  are  impeded  by  stiff 
londitudinal  springs  (see  last  chart  of  Fig.  14). 

The  longitudinal  springs  have  different  influence  on  the  first,  second  and  third  natural  frequencies 
of  vibration.  For  the  second  mode  of  vibration  large  crack  sliding  displacements  are  observed  in  the 
absence  of  bridging  mechanisms,  for  a  =  0  and  /}  -  0  (Fig.  15  and  first  chart  of  Fig.  12). 
Consequently,  when  present,  the  longitudinal  springs  will  be  significantly  engaged  and  their  influence 
on  the  frequency  ratio  will  be  strong  (Fig.  6).  For  the  third  mode  of  vibration,  on  the  other  hand, 
negligible  crack  sliding  displacements  are  observed  in  the  absence  of  crack  bridging  mechanisms  (Fig. 
15  and  first  chart  of  Fig.  13)  and  consequently  the  longitudinal  springs  will  not  have  a  strong  influence 
on  the  frequency  ratios  unless  the  stiffness  of  the  transversal  springs  is  simultaneously  increased  to 
prevent  the  opening  (Fig.  7).  In  this  case  the  beam  then  behaves  as  explained  above  for  the  constrained 
mode  case. 


9.  APPLICATION  TO  A  STITCHED  QUASI-ISOTROPIC  LAMINATE 

The  model  formulated  in  the  previous  sections  has  been  applied  to  investigate  if  a  through-thickness 
reinforcement  could  contrast  the  degradation  of  the  dynamic  properties  of  delaminated  beams.  To  this 
purpose,  a  48-ply  carbon/epoxy  quasi-isotropic  laminate,  of  thickness  7.2  mm,  stitched  with  glass  fiber 
tows  and  tested  in  [35-36]  has  been  considered.  The  laminate  had  a  reduced  longitudinal  Young’s 
modulus  of  49  GPa  and  a  shear  modulus  of  2.6  GPa.  The  stitches  were  spaced  according  to  a  square 
array  of  side  3.2  mm,  with  a  resulting  stitch  area  fraction  of  0.062. 

Experimental  shear  tests  have  been  performed  in  [35]  in  order  to  define  the  constitutive 
relationship,  shear  load  versus  crack  sliding  displacement,  of  a  single  stitch.  The  tests  were  carried  out 
loading  in  shear  small  notched  cubes  containing  just  one  or  two  stitches.  Experimental  curves 
representing  the  shear  force  applied  to  the  stitch  versus  the  relative  sliding  displacement  are  shown  in 
Fig.  16.  The  experimental  data  show  a  linearly  increasing  branch  followed  by  a  nonlinearly  decreasing 
one  at  a  large  value  of  the  sliding  displacement  corresponding  to  the  failure  of  the  stitch.  The  results 
presented  in  Fig.  16  have  been  used  to  define  the  constitutive  law  of  the  longitudinal  springs  of  the 
theoretical  model  by  assuming  an  equivalent  uniform  distribution  of  reinforcements  (see  right  ordinate 
axis  of  the  diagram  of  Fig.  16).  The  approximation  of  discrete  stitches  as  uniformly  distributed  springs 
is  correct  provided  the  stitch  density  is  high  enough  and  the  length  of  the  delamination  is  sufficiently 
long.  For  relatively  small  values  of  the  sliding  displacement  the  experimental  data  can  then  be 
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approximated  by  a  linear  bridging  law  describing  the  first  branch  of  the  experimental  curve.  The 
equivalent  normalized  stiffness  of  the  uniformly  distributed  springs  of  the  theoretical  model  is  /?=  0.01. 

Cuboidal  specimens  similar  to  those  tested  in  shear  were  also  tested  in  axial  tension  in  [36]  and  the 
stitches  showed  an  axial  stiffness  much  greater  than  the  shear  stiffness.  This  implies  that  the  value  of 
the  normalized  stiffness  of  the  transversal  springs,  a ,  is  much  higher  than  all  values  of  ac  presented  in 
Section  7  so  that  the  constrained  model  assumption  is  always  satisfied. 

Circles  have  been  used  to  highlight  the  solutions  corresponding  to  /?  =  0.01  in  Figs.  2-4.  The  data 
in  Fig.  2  show  that  for  a  mid-plane  delamination,  t2/f  =  0.5,  the  reinforcements  are  able  to  raise  the 
/ //„  ratio  from  0.926  to  0.992.  Figure  3  shows  that,  for  the  particular  geometry  considered,  P  =  0.01 
ensures  that  no  reduction  is  observed  in  the  first  natural  frequency  regardless  of  the  location  of  the 
delamination  along  the  beam.  Finally,  the  chart  of  Fig.  4  indicates  that  the  third  natural  frequency 
mostly  benefits  of  the  reinforcements:  the  ///0  ratio  is  brought  from  84.4  %  to  94.7%  of  its  intact 
value.  The  recovery  is  even  greater  for  longer  delaminations.  Results  corresponding  to  the  material 
tested  in  [35]  can  be  observed  in  the  diagrams  of  Figs.  5-10  for  P  =  0.01  and  a  corresponding  to  the 
constrained  model  (c.m.). 


10.  CONCLUSIONS 

A  theoretical  model  has  been  formulated  to  investigate  how  linear  elastic  bridging  mechanisms  can 
reduce  the  degrading  effects  of  delaminations  on  the  dynamic  response  of  beam-type  structures.  The 
model  is  based  on  the  first  order  shear  deformation  theory  for  laminated  plates  and  includes  the  effect 
of  the  rotational  inertia.  The  bridging  mechanisms  are  represented  by  uniformly  distributed  springs  that 
oppose  transversal  and  longitudinal  relative  displacements  between  the  surfaces  of  the  delamination. 

Assuming  no  extensional/bending  coupling  in  the  beam  constitutive  equations,  negligible  shear 
deformations  and  the  absence  of  opening  modes,  the  model  has  been  applied  to  analyze  cantilever 
beams  with  rectangular  cross  sections  and  delaminations  differently  located  in  the  longitudinal  and 
transversal  directions.  The  effectiveness  of  the  longitudinal  springs  on  reducing  the  natural  frequency 
shifts  has  been  demonstrated  through  die  transition  from  damaged  to  intact  beam  values  upon 
increasing  the  stiffness  of  the  springs. 

The  model  has  then  been  applied  to  a  geometrical  configuration  of  the  delaminated  beam  in  which, 
in  the  absence  of  bridging  mechanisms,  opening  modes  would  occur.  The  influence  of  both 
longitudinal  and  transversal  spring  stillness  on  the  first  natural  frequencies  and  corresponding  modal 
shapes  has  been  assessed.  It  has  been  shown  that,  by  increasing  the  transversal  longitudinal  spring 
stiffness  from  zero  to  relatively  high  values,  the  modal  shape  of  the  delaminated  beam  undergoes  a 
transition  during  which  the  opening  mode  progressively  disappears.  Simultaneously,  the  global  mode 
of  vibration  of  the  beam  changes:  from  a  first  to  a  second  mode  of  vibration  for  the  second  modal 
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shape  and  from  a  second  to  a  third  mode  of  vibration  for  the  third  modal  shape. 

The  model  has  been  formulated  for  composite  laminates  characterized  by  extensional/bending 
constitutive  coupling.  Results  have  been  presented  only  for  a  homogeneous  material  and  assuming 
negligible  shear  deformations.  It  is  expected  that  the  quantitative  results  will  be  modified  if  the  lay-up 
of  the  laminate  is  not  symmetric  about  the  mid-planes  of  the  beams.  However,  the  general  conclusions 
drawn  on  the  dynamic  response  of  delaminated  beams,  the  transitions  observed  on  varying  the  stiffness 
of  the  springs  and  the  influence  of  a  through-thickness  reinforcement  are  expected  to  hold  true. 

The  application  of  the  model  to  a  laminate  composite  reinforced  through  the  thickness  by  stitching 
and  previously  tested  in  [35],  has  shown  that  low  percentages  of  through-thickness  reinforcements  can 
substantially  improve  the  dynamic  response  of  delaminated  beams  and  restore  the  conditions  of  an 
intact  beam.  In  particular,  low  area  fractions  of  through-thickness  reinforcement  can  ensure  the  absence 
of  opening  modes  during  vibration  also  for  geometries  characterized  by  very  long  delamination  cracks 
close  to  the  surface  of  the  beam. 
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APPENDIX  I:  NUMERICAL  SOLUTION  OF  THE  BOUNDARY  EIGENVALUE  PROBLEM 

For  the  numerical  solution  of  the  boundary  eigenvalue  problem  formulated  in  Section  4  the  governing 
equations  (10-18)  are  first  transformed  into  an  equivalent  sets  of  first  order  ordinary  differential 
equations.  To  this  purpose,  the  following  notation  is  introduced: 
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where  it  is  understood  that: 


(Al.l) 


(A1.2) 


The  unknowns  quantities  of  the  problem  are  the  32  displacements  introduced  in  (Al.l,  A1.2)  and 
die  eigenvalue  <// .  The  system  is  completed  by  the  boundary  and  continuity  conditions  (19-24)  and  by 
artificially  adding  the  equation 


^»f  =  0  (A1.3) 

9 

which  merely  states  that  the  eigenvalue  is  independent  from  the  dimensionless  abscissa  £  [37]. 

The  system  of  governing  equations  and  related  boundary  and  continuity  conditions  define  a  two- 
point  boundary  value  problem  with  internal  boundary  conditions  which  are  solved  numerically  using  a 
relaxation  technique  [34,37].  The  dependent  variables  are  arbitrarily  numbered  in  order  to  define  the 
unknown  vector  y  =  (yj  >’2  " 1  }•  The  first  order  ordinary  differential  equations  (ODEs)  are  then 

replaced  by  approximate  finite  difference  equations  (FDEs)  on  a  grid  of  points  that  spans  the  interval  of 
interest.  The  present  problem  involves  N=  33  coupled  first-order  ODEs  approximated  by  just  as  many 
FDEs  on  a  mesh  of  M  points,  leading  to  a  total  of  Nx  M  unknowns. 

The  relaxation  method  provides  the  solution  by  starting  with  an  initial  guess  that  is  improved 
iteratively.  As  for  the  iteration  scheme,  the  multidimensional  Newton’s  method  provides  a  simple 
algebraic  matrix  equation  to  be  solved  for  the  N  corrections  Ay,-  at  the  M  grid  points.  Following  the 
procedure  described  in  [37]  for  two-boundary  value  problems  and  adapting  it  to  deal  with  the  presence 
of  the  internal  conditions,  the  matrix  of  the  algebraic  system  takes  a  special  “block  diagonal”  form  that 
allows  efficient  inversion.  More  details  on  the  numerical  solution  can  be  found  in  [33]. 


APPENDIX  II:  DERIVATION  OF  THE  GENERAL  SOLUTION 


The  general  integrals  (32-35)  of  Eqs.  (25a, b)  and  (26a, b)  are  straightforward.  The  following  parameters 
have  been  introduced  in  (32-35)  with  i  =  1  and  4: 


Si=4v, 


(A2.1) 
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In  order  to  determine  Eqs.  (36,37),  Eq.  (30)  has  to  be  solved  with  respect  to  wr,5 .  Substituting  the 
result  into  Eq.  (29),  after  differentiation  with  respect  to  £,  provides  the  following  ordinary  differential 
equation  in  the  unknown  w : 


(A2.2) 


where  the  quantities^,  8  and  </>  are  functions  of  the  eigenvalue  y/  and  the  stiffness  parameter  p  : 

X  =  2(y/-S p\  8  =  y/{y/-\6p~4%\  ^  =  48^(4 p~yi)  .  (A2.3) 

The  general  solution  of  Eq.  (A2.2)  corresponds  to  Eq.  (37),  where  the  quantities  al  are  the  solutions  of 
the  cubic  characteristic  equation  associated  to  the  differential  equation  (A2.2). 

Differentiating  Eq.  (30)  with  respect  to  £  and  substituting  the  result  into  Eq.  (29)  provides: 


«r=7 E-*»wk+M«s+m&,?]  > 


(A2.4) 


where  r ,  jj  and  v  are  given  by: 

r  =  24p(y/-4p\  n  =  Yip-yf,  v  =  48 {y/~p)  . 


(A2.5) 


Substituting  Eq.  (37)  into  Eq.  (A2.5)  provides  the  general  integral  (36),  where  the  quantities  f,  and  gt 
take  the  following  form: 


J!-a 

¥ 


V+M_48! t+2/>\  e,=-M 2L-M.-igr.-2fi 

kt  t  t  J  y/\t  r  t  j 


(A2.6) 


Finally,  rewriting  Eq.  (27)  as: 


u2  +y/u2  —  —2p\ 


\ 

“'+^J  ’ 


(A2.7) 


and  substituting  Eqs.  (36,37),  leads  to  the  general  integral  (38),  with  pt  and  qt  given  by: 


Pi  =  -2 q,=-2/7^  +  |j.  . 


(A2.8) 
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FIGURE  CAPTIONS 

Figure  1:  a)  Schematic  of  the  four-beam  model  for  a  delaminated  beam;  b)  Stress  resultants  and 
displacements. 

Figure  2:  Normalized  first  natural  frequency  of  vibration  as  a  function  of  the  normalized  stiffness  p  for 
different  transversal  locations  of  the  delamination  (constrained  model). 

Figure  3:  Normalized  first  natural  frequency  of  vibration  as  a  function  of  the  normalized  stiffness  P  for 
different  longitudinal  locations  of  the  delamination  (constrained  model). 

Figure  4:  Normalized  first,  second  and  third  natural  frequencies  of  vibration  as  a  function  of  the 
normalized  stiffness  P  for  a  mid-plane  delamination  centrally  located  (constrained  model). 

Figure  5:  Dependance  of  the  frequency  ratio  on  the  normalized  stiffness  P  for  the  first  natural 
frequency  (s.m.:  shear-bridged  model;  c.m.:  constrained  model). 

Figure  6:  Dependance  of  the  frequency  ratio  on  normalized  stiffness  P  for  the  second  natural 
frequency  (s.m.:  shear-bridged  model;  c.m.:  constrained  model). 

Figure  7:  Dependance  of  the  frequency  ratio  on  the  normalized  stiffness  p  for  the  third  natural 
frequency  (s.m.:  shear-bridged  model;  c.m.:  constrained  model). 

Figure  8:  Dependance  of  the  frequency  ratio  on  the  normalized  stiffness  a  for  the  first  natural 
frequency  (s.m.:  shear-bridged  model;  c.m.;  constrained  model). 

Figure  9:  Dependance  of  the  frequency  ratio  on  the  normalized  stiffness  a  for  the  second  natural 
frequency  (s.m.:  shear-bridged  model;  c.m.:  constrained  model). 

Figure  10:  Dependance  of  the  frequency  ratio  on  the  normalized  stiffness  a  for  the  third  natural 
frequency  (s.m.:  shear-bridged  model;  c.m.:  constrained  model). 

Figure  11:  Evolution  of  the  first  modal  shape  versus  the  normalized  stiffness  a ,  assuming  p~  0. 

Figure  12:  Evolution  of  the  second  modal  shape  versus  the  normalized  stiffness  a ,  assuming  p-  0. 

Figure  13:  Evolution  of  the  third  modal  shape  versus  the  normalized  stiffness  a ,  assuming  p=  0. 

Figure  14:  Evolution  of  the  third  modal  shape  versus  the  normalized  stiffness  P ,  assuming  a  = 
3.5*  1(T8. 

Figure  15:  Dependence  of  the  parameter  X  on  the  transversal  spring  stiffness  a ,  assuming  /?=  0. 

Figure  16:  Shear  load  versus  sliding  displacement  from  cuboidal  specimens  of  a  stitched  quasi¬ 
isotropic  laminate  containing  one  S-2  glass  fiber  tow  (after  [35]). 


23 


FIGURE  1 


L/t=50  -  Ldel/L  =  0,50  - 


FIGURE  3 


t2/t  =  0.5 


1.E-02  1.E-01  1.E+00 


Longitudinal  Spring  Stiffness  -  Beta 


FIGURE  5 


Longitudinal  Spring  Stiffness  •  Beta 


FIGURE  6 


FIGURE  7 


/ 


FIGURE  9 


Transversal  Spring  Stiffness  -  Alpha 


FIGURE  10 


FIGURE  13 


oo  co  Tt-  cm  o  ca 


B  =  10'4  C'fO  =  0.800 


FIGURE  14 


FIGURE  15 


60 


Sliding  displacement  (mm) 


FIGURE  16 


Equivalent  bridging  traction 
itb  0.099  stitches/mm2  (MPa) 


reinforcement  percentage  and  fiber  orientation  for  the  dif¬ 
ferent  layers.  This  type  of  structure  should  allow  for  a 
significant  weight  reduction  in  comparison  with  similar 
metallic  shields,  while  the  overall  energy  absorption  is  ex¬ 
pected  to  satisfy  the  application  requirements  that,  t»Wng 
into  account  the  train  mass,  are  very  severe. 

The  first  phase  of  die  research  activity  has  been  devoted 
to  the  experimental  characterization  of  the  proposed  mate¬ 
rial,  with  particular  attention  to  the  interface  between  the 
composite  layer  and  the  foam.  To  this  aim  a  number  of 

quasi-static  and  impact  tests  have  been  run  both  on  tensile 

and  bending  specimens  of  the  composite  laminates  and  on 
square  plates  of  the  whole  sandwich. 

The  second  phase  of  toe  research  activity  has  been  de¬ 
voted  to  toe  finite  element  modeling  of  toe  characteriza¬ 
tion  tests  run  in  the  first  phase.  The  purpose  of  this  second 
phase  was  to  assess  toe  capability  of  toe  FE  code  mate¬ 
rial  model  to  simulate  toe  experimental  tests  and  to  define 
toe  input  data  for  toe  material  constants  required  by  toe 
chosen  material  model  The  third  and  final  phase  of  toe 
research  activity  has  been  dedicated  to  the  study  of  the 
shield  structure  impact  behavior  by  means  of  finite  ele¬ 
ment  analysis.  The  FEM  simulation  was  performed  with 
toe  explicit  code  PAM-Crash. 
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The  paper  deals  with  different  aspects  of  the  mechan- 
ical  behavior  of  delaminated  composites  in  toe  pres¬ 
ence  of  bridging  mechanisms  such  as  those  developed  by 
a  through-thickness  reinforcement  (stitching,  ^-pinning 
or  weaving).  We  will  present  our  recent  work  on  the 
problems  of  mode  I  and  mixed  mode  large  scale  bridg¬ 
ing  delamination  and  toe  dynamic  response  of  through- 
thickness  reinforced  delaminated  beams.  The  problem  of 
mode  I  and  mixed  mode  large  scale  bridging  delamina¬ 
tion  is  treated  using  toe  weight  function  method  in  terms 
of  stress  intensity  factors  at  toe  crack  tip.  The  model  ac¬ 
counts  for  toe  ortootropy  of  toe  material  and  allows  us  to 
investigate  problems  characterized  by  very  short  cracks 
or  bridging  mechanisms  acting  close  to  toe  crack  tip.  The 
influence  of  toe  through-thickness  reinforcement  on  toe 


fracture  behavior  is  investigated  and  unusual  phenomena 
of  crack  arrests  and  crack  propagation  in  the  presence  of 
regions  of  contact  that  have  been  observed  in  the  exper¬ 
iments  are  predicted  theoretically  and  explained  (1J.  A 
beam  theory  model  has  been  formulated  to  analyze  the 
role  played  by  linear  bridging  mechanisms  in  reducing  toe 
delammation-mduced  degradation  of  the  dynamic  proper¬ 
ties  of  laminated  beams.  The  model  follows  toe  formula¬ 
tion  proposed  in  [2]  for 

delaminated  beams  in  toe  absence  of  bridging  mecha¬ 
nisms.  An  application  of  toe  model  to  a  typical  stitched 
carbon  epoxy-laminate  shows  that  a  low  area  fraction  of 
through-thickness  reinforcement  can  substantially  reduce 
toe  degradation  of  toe  natural  frequencies  of  the  structure. 

[1]  Massabd,  R.,  Brandinelli,  L.,  and  Cox  B.N.,  (2002) 
Mode  I  Weight  Functions  for  an'  Orthotropic  Double  Can¬ 
tilever  Beam,  accepted  for  publication  in  Ini  J.  Engineer¬ 
ing  Science. 

[2]  Luo  H.,  Hanagud  S.  (2000),Dynamics  of  delaminated 
beams,  Intcmatinal  Journal  of  Solids  and  Structures,  Vol 
37,  pp.1501-1519. 

Hie  authors  were  supported  by  the  European  Research  Office  of  die  US 
Army,  Contract  #  N68171-01-M-5909.  and  by  the  Itafian  Department 
for  the  University  and  for  Scientific  and  Ibdmological  Research. 


Contact  Between  a  Cylindrical 
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In  this  paper  we  analyze  toe  problem  of  indentation  of 
a  sandwich  beam  by  a  rigid  cylindrical  indentor.  The  core 
is  assumed  to  be  functionally  graded  with  stiffness  prop¬ 
erties  varying  through  toe  thickness.  The  contact  prob¬ 
lem  is  solved  using  toe  assumed  contact  stress  distribution 
method.  Relations  are  derived  for  load-contact  length, 
load-indentation,  and  maximum  stresses  and  strain^  in  thft 
core  for  functionally  graded  cores  as  well  as  homoge¬ 
neous  core.  The  results  will  be  useful  In  answering  toe 
question  if  FG  cores  have  an  advantage  over  homoge¬ 
neous  cores  when  a  sandwich  panel  is  subjected-  to  low- 
velocity  impact 

The  core/  face  sheet  delamination  is  a  major  concern  in 
sandwich  laminate  construction.  The  stiffness  discontinu¬ 
ity  at  toe  face  sheet  and  core  interface  results  in  a  large 
increase  in  shear  stresses.  While  toe  core  material  itself 


15th  ASCE  Engineering  Mechanics  Conference 
June  2-5, 2002,  Columbia  University,  New  York,  NY 


FREE  VIBRATIONS  OF  THROUGH-THICKNESS  REINFORCED 
DELAMINATED  BEAMS 
Luigi  Brandinelli 1  and  Roberta  Massabo 2 

ABSTRACT 

The  dynamic  behavior  of  laminated  beam-type  structures  presenting  delaminations  is  characterized  by 
stiffness  degradation,  natural  frequency  shifts  and  opening  modes  during  vibration.  This  paper  investigates 
the  influence  of  a  through-thickness  reinforcement  on  the  natural  frequencies  of  through-width  delaminated 
beam.  A  model  based  on  the  theory  of  bending  of  beams  is  formulated  where  the  reinforcements  are 
modeled  as  a  uniform  distribution  of  linear  springs  able  to  react  to  longitudinal  and  transversal  relative 
displacements  between  the  surfaces  of  the  delamination.  The  natural  frequencies  of  a  delaminated 
cantilever  beam  are  predicted  for  different  lengths  of  the  delamination.  The  effectiveness  of  the  through 
thickness  reinforcement  on  the  natural  frequencies  is  demonstrated  through  the  transition  from  damaged 
(unbridged)  to  intact  beam  values  upon  increasing  the  stiffness  of  the  springs.  An  application  of  the  model 
to  a  typical  stitched  carbon-epoxy  laminate  for  aeronautical  applications  shows  that  low  percentages  of 
through-thickness  reinforcements  can  substantially  improve  the  dynamic  response  of  delaminated  beams. 

Keywords:  Through-thickness  reinforcement,  delamination,  natural  frequencies,  beam  theory, 
dynamics. 

INTRODUCTION 

Delamination  is  one  of  the  most  common  failure  mechanisms  of  laminate  composite 
materials.  Besides  being  known  as  major  cause  of  stiffness  and  strength  degradation, 
delamination  is  also  known  to  cause  a  change  in  the  dynamic  behavior.  Extensive  experimental 
and  theoretical  investigations  have  demonstrated  that  delaminations  can  cause  shifts  and/or  extra 
natural  frequencies  which  can  approach  the  exciting  frequencies.  Additionally,  delaminations  can 
cause  changes  in  modal  shapes.  Conditions  causing  the  opening  of  delaminated  surfaces  in  modal 
shapes,  so-called  “opening  modes”,  have  been  experimentally  detected  and  theoretically 
reproduced  (Wang  et  al.,  1982;  Shen  and  Grady,  1992;  Luo  and  Hanagud,  2000). 

The  efficacy  of  the  through-thickness  reinforcement  technology  (stitching,  z-pinning, 
weaving)  in  enhancing  the  delamination  strength  and  damage  tolerance  of  composite  laminates 
has  been  abundantly  documented  in  the  last  twenty  years.  In  particular,  it  has  been  proven  how 
small  percentages  of  reinforcements  can  increase  interlaminar  fracture  toughness,  shear  and 
compressive  strength  after  impact  (Dransfield  et  al.,  1994).  On  the  other  hand,  to  the  knowledge 
of  the  authors,  the  efficacy  of  through-thickness  reinforcements  in  limiting  the  delamination- 
induced  degradation  of  the  dynamic  properties  has  not  yet  been  investigated.  It  is  reasonable  to 
assume  that  through-thickness  reinforcements  can  reduce  the  degradation  of  the  dynamic 
properties  by  opposing  the  relative  displacements  between  the  delamination  surfaces. 


1  Dept,  of  Structural  end  Geotechnical  Engineering,  University  of  Genoa,  Italy.  E-mail:  brandinelli@diseg.unige.it 
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A  simple  model  based  on  laminated  beam  theory  has  been  devised  in  order  to  study  the  role 
played  by  through-thickness  reinforcements  on  the  dynamic  response  of  delaminated  composite 
beam.  Following  the  approach  described  in  Luo  and  Hanagud  (2000),  a  through-width 
delaminated  beam  has  been  modeled  by  using  four  elementary  beams  connected  at  the 
delamination  edges.  The  action  of  the  reinforcement  has  been  represented  by  a  uniform 
distribution  of  linear  elastic  springs  able  to  oppose  longitudinal  and  transversal  relative  crack 
displacements  (Massabd  and  Cox,  1999).  The  constitutive  law  of  the  springs  describes  the  load 
transfer  mechanisms  from  the  reinforcements  to  the  surrounding  material  and  is  assumed  as  a 
model  input.  Although  a  more  general  formulation  of  the  model  could  be  based  on  a  nonlinear 
spring  constitutive  law,  in  the  present  work  a  linear  relationship  has  been  assumed,  which  is  an 
acceptable  assumption  for  most  through-thickness  reinforced  laminates  provided  the  crack 
displacements  are  not  too  high  (Cox  et  al.,  1997).  The  model  assumes  bending/extensional 
coupling  in  the  laminate  constitutive  description,  and  accounts  for  the  effect  of  inertia  in  the 
moment  dynamic  equilibrium  equation  of  the  beam.  Shearing  deformation  effect  is  neglected. 
Through  this  approach,  delaminations  of  different  lengths  and  at  different  sites,  both  thickness- 
wise  and  length-wise,  can  be  analyzed. 

The  model  has  then  been  used  to  predict  the  effect  of  through-thickness  reinforcements  for  a 
wide  range  of  geometric  configurations,  reinforcement  stiffness  values  and  material  parameters 
of  the  composite  laminate.  Finally,  an  application  of  the  model  is  presented  to  simulate  the 
dynamic  response  of  a  stitched  carbon-epoxy  laminate  for  aeronautical  applications.  In  this  case, 
the  appropriate  stiffness  value  for  the  linear  elastic  springs  has  been  deduced  from  the 
experimental  data  in  Cox  et  al.  (1997). 

THEORETICAL  MODEL 

The  model  refers  to  the  scheme  shown  in  Fig.l  representing  a  laminated  beam  of  length  L 
and  thickness  t  with  a  through-width  delamination  of  length  Ldc\.  Linear  bridging  mechanisms 
acting  along  the  surfaces  of  the  delamination  are  represented  in  the  model  by  means  of  uniformly 
distributed  linear  elastic  springs,  with  stiffness  £v  in  the  transverse  direction  and  kx  in  the 
longitudinal  direction.  Plane  strain  conditions  parallel  to  plane  x~y  are  assumed  and  the  problem 
is  solved  by  applying  the  theory  of  bending  of  beams. 

The  delamination  identifies  four  beams,  for  /  =  1,  ...4,  of  length  Lit  thickness  fi}  distance  of 
the  centroid  from  the  delaminated  surfaces  h\  and  cross  sectional  area  and  moment  of  inertia  S\ 
and  I,.  The  generalized  displacements  are  the  transverse  displacement,  Wj,  the  axial  displacement, 
Mi,  and  the  bending  rotation,  $,  and  the  stress  resultants  per  unit  width  are  the  normal  force,  N-ti 
shear  force,  Q,  and  bending  moment,  Mh  for  i  =  1,..,  4  (Fig.  2). 

Equations  of  motion  and  constitutive  equations  are  defined  for  the  four  beams  and  a 
boundary  value  problem  is  set  up  in  terms  of  the  generalized  displacements  (see  Appendix  I). 
The  problem  generally  requires  a  numerical  solution.  However,  a  closed-form  solution  can  be 
obtained  based  on  the  following  assumptions: 

1)  mid-plane  delamination; 

2)  homogeneous  material; 

3)  absence  of  opening  modes:  the  transverse  displacements  of  the  two  delaminated  beams 

coincide,  that  is  w2  =  w3  =  w . 

Assumption  2  is  evidently  satisfied  not  only  by  all  isotropic  and  orthotropic  homogeneous 
materials  but  also  by  unidirectional  lay-ups  made  of  the  same  type  of  orthotropic  lamina  and 
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laminates  with  many  layers  and  no  bending/extensional  coupling.  Quasi-isotropic  laminates 
typically  used  for  aeronautical  applications  and  cross-plies  [0/90]  belong  to  this  category. 
Assumption  3,  is  usually  satisfied  if  the  delamination  is  not  too  close  to  the  surface  of  the  beam, 
as  observed  by  Luo  and  Hanagud  (2000),  and  is  commonly  referred  to  as  “constrained  model” 
assumption. 


FIG.  1:  Scheme  of  the  four-beam  model  for  a  delaminated  beam. 


FIG.  2:  Sign  convention  for  stress  resultants  and  displacements. 


Assuming  a  free  harmonic  vibration  characterized  by  a  single  angular  frequency  ©, 
Wj  (x,  t)  =  u{  (x)eJO*  and  (x, /)  =  wt  (x)ej0*  for  i  =  1,  ..4,  the  problem  is  defined  in  terms  of  the 
relative  axial  displacement  amplitude,  ur  =«2  -u3,  and  the  transverse  displacement  amplitude, 
w,  of  the  beams  2  and  3  and  by  the  axial  and  transverse  displacement  amplitudes  of  the  beams  1 
and  4. 

It  is  convenient  to  formulate  the  problem  in  a  dimensionless  form  introducing  the  normalized 
coordinate  ^  =  x  / /  and  the  normalized  amplitude  displacements  ur  =  ur  It  and  w  =  wit ,  for 
the  delaminated  beams,  and  u.  =  ujt  and  w.  =  wt  It ,  with  i  =  1,4,  for  beams  1  and  4.  The  other 
dimensionless  variables  are  the  spring  stiffness  p  and  the  eigenvalue  \jf  defined  as: 

2  2 

P=i k,tlEx  and  y/  = P^-  —  (1) 

E\ 

where  E{  is  the  effective  longitudinal  Young’s  modulus  of  the  beam  and  p  is  the  density. 

The  boundary  eigenvalue  problem  is  then  defined  by  the  following  equations: 
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i  =  1...4 


|«i»K+V^i  =0 

‘U„^+(H'~4P)ut+2Pw,i=0  (2) 

“(1 2j 8  -  y/)w,%  -48vav  =  24  pit, ,  ? 

The  continuity  conditions  at  the  first  delamination  edge  £  =  L^jt  are: 

#i  =  *.  w„f  =  w)f,  u2  =Ul+-^-Wi,(, 

Mx  =  M2+M3  -(-N2  +  *-N3,  Ql=Q2+Q3, 

4  4 

while  those  at  the  other  edge  %=  L2/t  (with  L2=Ll+LJel)  are: 

w4=w,  w4,(  =  w,(,  u2=u4+-^w4,(, 

M4  =  M2  +m3  ~n2  +^N3,Q4=Q2+Q3, 

The  boundary  conditions  of  a  cantilever  beam  are: 

Wj=0,  #^  =  0  at£=0 

M4=  0,  &  =  0  at  |=1. 

The  solution  procedure  is  outlined  in  Brandinelli  and  Massab6  (2001).  The  natural 
frequencies  and  the  corresponding  modal  shapes  of  the  delaminated  beam  are  calculated  through 
the  analytical  integration  of  the  eigenvalue  problem  (2).  It  is  worth  noting  that  for  very  low  and 
very  large  values  of  /?,  the  normalized  stiffness  of  the  ligaments,  a  high  numerical  precision  is 
required  in  the  solution. 

Under  these  assumptions,  the  proposed  model  has  been  validated  by  referring  to  the 
asymptotic  frequency  values  obtained  by  letting  the  spring  stiffness  vanish  (thus  simulating  the 
absence  of  reinforcement)  and  indefinitely  grow  (absence  of  delamination).  For  a  vanishing 
spring  stiffness,  the  experimental  tests  of  Shen  and  Grady  (1992),  on  delaminated  cantilever 
beams  without  a  through-thickness  reinforcement,  have  been  successfully  reproduced 
(Brandinelli  and  Massabd,  2001).  Luo  and  Hanagud  (2000)  used  the  same  tests  to  validate  their 
beam  model.  For  very  high  spring  stiffness,  the  frequency  values  predicted  by  the  model 
approach  the  values  corresponding  to  an  intact  beam. 


A  A  »  A 

«3  =M‘  4LWi,< 
Nt=N2  +  N3 

"3  =u4—^-w4ti 

n4=n2+n3 
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If  assumption  1  (mid-plane  delamination)  is  removed,  the  procedure  outlined  above  is  no 
longer  applicable.  It  is  still  possible  however  to  proceed  with  a  numerical  integration  of  the 
differential  eigenvalue  problem  described  in  Appendix  I,  based  on  relaxation  techniques  (Keller, 
1968).  The  system  of  differential  equations  (A4)  is  approximated  by  finite-difference  equations 
and  the  multidimensional  Newton-Raphson  method  is  used  for  the  iteration  scheme. 


FREE  VIBRATIONS  OF  THROUGH-THICKNESS  REINFORCED  BEAMS 

The  model  described  in  the  previous  section  has  then  been  used  to  systematically  determine 
the  role  played  by  through-thickness  reinforcements  in  reducing  the  delamination-induced 
degradation  of  the  dynamic  properties  of  through-width  delaminated  cantilever  beams. 

The  dimensionless  formulation  of  the  problem  allows  the  definition  of  the  calculated  natural 
frequency/  of  the  delaminated  beam,  normalized  with  respect  to  the  corresponding  intact  beam 
value/,  as  a  junction  of  the  dimensionless  stiffness  of  the  ligaments  /3  —ktt!Ex,  the 

normalized  length  of  the  delamination,  LdeljL,  the  normalized  length  of  the  beam,  Ljt  and  a 
parameter  describing  the  position  of  the  delamination  along  the  length,  {Lx  +  Ldel  j 2)/ L .  The 
chart  shown  in  Fig.3  refers  to  the  first  natural  frequency,  a  mid-plane  delamination,  Ljt =  50 
and  (ij  +  Ldel[2)lL  =  0.5.  The  three  curves  in  the  diagram  refer  to  three  different  values  (0.25, 
0.50  and  0.75)  of  the  L^/Z,  ratio.  For  very  small  values  of  p  the  curves  tend  to  the  normalized 

frequencies  of  the  delaminated  beam  in  the  absence  of  through-thickness  reinforcement.  As 
expected,  the  curves  show  a  transition  from  the  damaged  and  unbridged  values  to  the  values 
corresponding  to  an  intact  beam  upon  increasing  the  elastic  spring  stiffness. 

The  diagram  also  shows  the  values  of  ///0  corresponding  to  a  typical  stitched 

carbon/epoxy  laminate  composite  for  aeronautical  applications  (Cox  et  al.,  1997).  The  composite 
has  a  quasi-isotropic  lay-up  (48-plies  [45/0/-45/90])  and  a  thickness  of  7.2  mm.  The  reduced 
longitudinal  Young’s  modulus  is  Ei  =  49  GPa,  and  in-plane  Poisson’s  ratios  are  V12  —  V21  =  0.3. 

I - - - -  The  laminate  was  stitched  with  doubled 

j  First  Natural  Frequency  -  Lte=50  3550  denier  S-2  glass  fiber  tows  on  a 

square  array  of  side  3.2  mm;  the 
individual  stitch  area  is  0.636  mm2  and 
the  total  stitch  area  fraction  is  cs  = 
0.062. 

Experimental  shear  tests  have  been 
performed  by  Cox  et  al.  (1997)  in  order 
to  define  the  constitutive  relationship, 
load  versus  crack  sliding  displacement, 
of  a  single  reinforcement.  The 
—  experimental  data  showed  a  linearly 

1.&O6  IE-05  i.E-04  1.&03  1.&Q2  l.&oi  i.&oo  increasing  branch  followed  by  a 

0843  nonlinearly  decreasing  one  at  large 

FIG.  3.  Normalized  first  natural  frequency  dlsP^c®ment  corresponding  to 

as  a  function  of  the  normalized  stiffness  of  the  Jllur\of  the  reinforcement.  Based 
...  ......  .  o  ,  .  ...  on  this  observation  the  experimental 

the  reinforcement,  /J  =  k,t/Et ,  for  three  m  cm  be  approximated  with  a  linear 

normalized  delamination  lengths.  bridging  law  describing  the  first  branch. 

The  equivalent  normalized  stiffness  for 


1.&03  1.&Q2 

Beta 


FIG.  3.  Normalized  first  natural  frequency 
as  a  function  of  the  normalized  stiffness  of 
the  reinforcement,  f3  =ktt/El ,  for  three 
normalized  delamination  lengths. 
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the  spring  distribution  of  the  model  has  been  calculated  as  being  equal  to  j3«  0.01. 

In  the  diagram,  a  circle  is  used  to  highlight  the  solutions  corresponding  to  this  case.  They 
show  that  low  area  fractions  of  through-thickness  reinforcements  can  substantially  reduce  the 
delamination-induced  degradation  of  the  dynamic  properties  of  the  beam.  For  a  delamination 
with  LdeljL  =  0.5,  the  through  thickness  reinforcement  takes  the  value  of  the  first  natural 

frequency  from  0.93  f0  back  to  0.994 f0  with  a  7%  relative  increment.  Moreover,  the  effectiveness 
of  the  through-thickness  reinforcements  on  the  frequency  heightens  upon  increasing  the 
normalized  length  of  the  delamination.  Analyses  carried  out  for  different  values  of  the 
normalized  length  of  the  beam  Ljt  show  that  the  effectiveness  of  the  reinforcement  increases 
upon  raising  the  Ljt  ratio. 


FIG.  4.  Normalized  second  and  third  natural  frequencies  as  a  function  of  the 
normalized  stiffness  of  the  reinforcement,  P=ktt/El9  for  three  different 
normalized  delamination  lengths. 

Figure  4  shows  diagrams  of  the  second  and  third  natural  frequencies.  As  already  noted  in 
experimental  tests  on  damaged  composites  (Kimpara  et  al.  1998)  the  damage  primarily  affects  the 
frequencies  of  high-order  vibration  modes.  For  the  case  L^jL  =  0.5,  for  instance,  the  damaged 

value  (unbridged)  of  the  third  natural  frequency  is  63%  of  the  corresponding  intact  value,  while 
the  damaged  value  of  the  first  frequency  is  93%  of  the  intact  value.  The  model  proposed  here  also 
shows  how  the  influence  of  the  bridging  mechanisms  acting  along  the  crack  is  more  evident  on 
the  frequencies  of  high-order  modes.  For  the  stitched-carbon  epoxy  laminate  described  above  and 
Ldet/L  =  0.5,  for  instance,  the  through-thickness  reinforcements  induce  a  43%  relative  increment 
on  the  damaged  value  of  the  third  frequency  and  only  a  7%  increment  on  the  first  frequency.  This 
observation  can  be  useful  for  the  design  of  dynamic  experimental  tests  to  investigate  the 
properties  of  through-thickness  reinforcements  through  the  solution  of  an  inverse  model. 

CONCLUSIONS 

A  theoretical  model  has  been  proposed  to  evaluate  how  the  translaminar  reinforcement  can 
limit  the  effects  of  delaminations  on  the  dynamic  response  of  beam  structures.  The  model  is 
based  on  laminated  beam  theory  and  includes  the  effect  of  the  inertia  in  the  beam’s  moment 
dynamic  equilibrium  equations.  Effect  of  shearing  deformation  is  neglected.  The  action  of 
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reinforcements  is  simulated  through  a  uniform  distribution  of  linear-elastic  springs  that  are  able 
to  counteract  the  relative  transversal  and  longitudinal  displacements  between  the  surfaces 
delineating  the  delamination.  Assuming  no  extensional/bending  coupling  in  the  beam  constitutive 
equations  and  absence  of  opening  modes,  the  model  has  been  applied  to  the  case  of  a  delaminated 
cantilever  beam  for  different  delamination  lengths.  The  effectiveness  of  the  through-thickness 
reinforcement  on  reducing  natural  frequency  shifts  has  been  demonstrated  through  the  transition 
from  damaged  to  intact  beam  values  upon  increasing  the  stiflhess  of  the  springs.  The  application 
of  the  model  to  a  laminate  composite  reinforced  through  the  thickness  by  stitching  and  previously 
tested  by  the  authors,  has  shown  that  low  percentages  of  through-thickness  reinforcements  can 
substantially  improve  the  dynamic  response  of  delaminated  beams. 
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APPENDIX  I. 

The  differential  equations  of  axial,  transversal  and  rotational  motion  for  beams  1  and  4  can  be 
simply  written  as 

N;,x  =  pS;U.„tt  Q^  =  ps^a  -Qi+M„x  =  pl^>i,u  i  =  l,4  (Al) 
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where  partial  differentiation  with  respect  to  the  spatial  variable  x  and  time  t  is  denoted  by  a 
comma.  For  the  delaminated  beams  the  equations  of  motion  include  the  contribution  of  the  action 
performed  by  the  springs: 

#  2  >x  +  *«  («3  ~  h3<P  3  “  «2  “  h2(p  2  )=  pS  2U  2 
N  J.X  -*«  («3  -  h3(p  3  ~  U2  -  h2<p  2  )=  pS  3«3  ,tt 
2  2  *x  +*  V  (W  3—  W2)=P^2W2»tt 

Q3>x~ky  (w3  -  w2)=  pS3w3,n  (A2) 

—  02  +M2,x+kth2(u3  ~  <p  3h3  —  U  2  ~  ?  2^2  )=  P I  2*P  2  »  tt 

—  Q  3  +  M3ix+kih3($i3  -  <p  3  A  3  -  m  2  —  ?  2^2  )=  P?  3^3*# 

The  constitutive  equations  of  the  laminated  beam  can  be  written  in  the  following  form: 

Nl=A»i>*+Bi<pix  (A3) 

where  4  and  A  represent  the  axial  and  bending  stiffness  respectively,  and  is  the  extension/ 
bending  coupling  coefficient.  Using  the  relation  (pi  =  -w,  x ,  and  introducing  the  quantities 

/Wj  =  pSj ,  Jj  =  p/j ,  the  governing  equations  of  the  model  are 

^tWi  >xxx“^i^i  ’tt 

4W2m-^W2.xxx^t(«3+V3.x-^+^2.x)=W2«2.tt 
-  4«3»xx-^3»xxx^t(«3+^»x-^+V2»x)=»h«3.tt 

A^2m-^“2»xxx^A("3»x^^.xx^2»x-^M’2.xx)-‘/2W2»xxtt-*v(^  ~W2)+m2W2m  =0  (A‘4) 

A^3.xx-AW3.xxx^t^(»3.x^W3>xx^2.x^W2.xx)-^3.xxtt-^v(w3  =0 
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